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1.  INTRODUCTION 


Strong  blast  waves  are  produced  in  gaseous  media  by  a  sudden  depo¬ 
sition  of  large  amounts  of  energy  in  a  relatively  small  region.  Such 
events  may  be  approximated  mathematically  by  a  point  source  wave.  If 
one  assumes  in  addition  an  ideal  gas  and  an  initial  pressure  negligible 
compared  to  the  shock  pressure,  then  the  flow  field  can  be  described  in 
closed  form  as  a  self-similar  solution  of  the  flow  equations^ » ^2 3 4  5 * 7 ^ .  A 
derivation  of  such  closed  form  solution,  including  spherical,  cylindrical 
and  planar  blast  waves,  is  given,  e.g.,  by  Sedov^  and  Rouse  .  The  so¬ 
lutions  have  been  supplemented  by  Laporte  and  Chang^5^  who  derived  closed 
form  expressions  also  for  the  particle  paths  and  for  the  Mach  lines. 


2 

H.A.  Bethe,  K.  Fuchs ,  J.  von  Neumann,  R.  Peierls,  and  W.G.  Penny ,  U.S. 
Atomic  Energy  Commission  Report  AECD-2860  (1944). 


2 

J.L.  Taylor,  " An  Exact  Solution  of  the  Spheirical  Blast  Wave  Problem", 
Phil.  Mag.,  Vol  46  (1955). 

3 

R.  Latter,  "Similarity  Solution  for  a  Spherical  Shock  Wave",  Journal  of 
Applied  Physics,  Vol  26,  No.  8  (1955). 


4 

L.I.  Sedov,  "Similarity  and  Dimensional  Methods  in  Mechanics" , 
Press,  New  York  (1959). 


Academic 


5 C.A .  Rouse,  "Theoretical  Analysis  of  the  Hydrodynamic  Flow  in  Exploding 
Wire  Phenomena,  in  " Exploding  Wires",  W.G.  Chace  and  K.M.  Howard,  edts.. 
Plenum  Press,  New  York,  pp.  227-263  (1959). 

Q 

0.  Laporte  and  T.S.  Chang ,  " Exact  Expressions  for  Curved  Characteristics 
Behind  Strong  Blast  Waves ”,  U.S.  Army  Ballistic  Research  Laboratory 
Contractor  Report ,  BRL-CR-SO  (January  1971).  (AD#722777) 


7 

0.  Laporte  and  T.S.  Chang ,  n Curved  Characteristics  Behind  Blast  Waves”, 
Physics  of  Fluids ,  Vol  15,  pp.  502-504  (1972) . 
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The  closed  form  solutions  are  valuable  approximations  of  physical 
phenomena,  and  they  have  been  used  as  test  cases  for  numerical  solvers 
of  flow  governing  equations*^  9.  The  evaluation  of  the  rather  cumbersome 
formulas  is,  however,  not  trivial  because  the  solution  is  expressed  in 
parameter  form.  Therefore,  also  several  tables  of  strong  blast  wave 
solutions  have  been  published^’^jlOjll.  The  use  of  such  tables  is,  of 
course,  limited  to  the  particular  set  of  parameters  chosen  by  the  authors 
of  the  tables.  Also,  because  tables  are  expressed  in  form  of  dimensionless 
variables,  they  must  be  transformed  before  application  to  suit  each  par¬ 
ticular  case.  Not  all  authors  of  the  published  tables  provide  sufficient 
and  clear  instructions  as  to  how  the  transformations  should  be  accomplished, 
nor  are  the  tables  of  various  authors  standardized.  Finally,  the  particle 
path  and  Mach  line  solutions  of  Laporte  and  Chang  have  not  been  tabulated 
at  all. 


In  order  to  make  the  closed  form  solutions  of  strong  blast  waves 
easily  available,  we  have  developed  a  package  of  computer  programs  that 
produces  the  solution  for  any  dimension  and  combination  of  parameters, 
within  limits.  The  present  report  is  a  description  of  the  programs  and 
of  their  use. 

In  Section  2  of  this  report  we  shall  summarize  the  formulas  needed 
for  flow  computation.  We  shall  not  give  a  derivation  of  the  formulas 
which  can  be  obtained,  e.g.,  from  References  4  and  6,  and  which  is  based 
on  the  following  assumptions: 

1.  the  flowing  medium  is  an  ideal  gas; 

2.  the  initial  pressure  in  the  ambient  gas  can  be  neglected; 


8P.C.  Chou  and  R.R.  Karpp ,  "Solution  of  Blast  Waves  by  the  Method  of 
Characteristics ",  Drexel  Institute  of  Technology  Report ,  DIT  Report 
No.  125-7,  (1965). 

9  J.A.  Schmitt,  "A  Finite  Element  Method  and  Corresponding  Pilot  Computer 
Code  for  Hyperbolic  Systems  of  Equations  in  Two  Spatial  Dimensions  and 
Time  Applied  to  Unsteady  Gas  Flows",  U.S.  Army  Ballistic  Research 
Laboratory  Report,  BRL-R-2017  (September  1977) . (AD#A045703) 

10N.  Gerber  and  J.M.  Bartos,  "Tables  of  Cylindrical  Blast  Functions  for 
y  =  5/3  and  y  =  7/5",  U.S.  Army  Ballistic  Research  Laboratory  Report, 
BRL-R-1376  (October  1961 ) . (AD#663821) 

Hj.W.  Goresh  and  R.G.  Dunn,  "Tables  of  Blast  Wave  Parameters,  I  Spherical 
Explosions",  Aerospace  Research  Laboratory  Report,  ARL-69-0011  (January 
1969). 
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3.  the  energy  is  released  instantaneously  at  a  point,  in  a  line  or 
in  a  plane. 

For  applications  of  the  results  to  real  life  situations,  the 
assumptions  mean  that  the  formulas  may  describe  accurately  events  that 
are  neither  too  far  from  the  explosion  (assumption  2) ,  nor  too  close  to 
it  (assumption  3) . 

In  Section  3  we  shall  describe  the  use  of  the  computer  programs 
for  strong  blast.  Some  examples  of  calculations  will  be  presented  in 
Section  4. 

In  Section  5  we  derive  quantitative  conditions  under  which  the 
strong  blast  solutions  can  provide  approximate  information  about  pressures 
generated  by  real  life  explosions  in  air. 

2.  THEORETICAL  BACKGROUND 


2.1.  Shock  Formulas 


Let  xs  be  the  location  of  the  shock,  i.e.,  the  distance  of  the  shock 

from  the  explosion,  and  let  t  be  the  time  elapsed  after  the  explosion. 

Using  dimensional  arguments^ or  group-theoretical  discussions^  one  can 

derive  the  following  relation  between  x  and  t: 

s 

/  E  \l/ (2+n)  2/ (2+n) 

xs  =  t  .  (2.1) 


In  Eq.  (2.1),  K(n,y)  is  a  proportionality  factor  that  depends  on  the 
ratio  of  specific  heats,  y,  and  on  the  dimension  n  of  the  event  (n  =  1 , 

2 ,  or  3  for  planar,  cylindrical,  or  spherical  symmetry,  respectively), 

Po  is  the  density  of  the  ambient  gas,  and  E  is  the  energy  released,  per 
unit  area  for  a  planar  wave,  per  unit  length  for  a  cylindrical  wave,  and 
the  total  energy  released  for  a  spherical  wave. 

The  proportionality  factor  K(n,y)  is  determined  from  the  requirement 
that  the  total  energy  contained  in  the  blast  wave  must  be  equal  to  the 
total  energy  released,  i.e.. 


E 

o 


n/ 2 

TT 


X  (t) 

/s  v  J 

p(x,t)  ^U2(x,t)  + 


p(x,t)\ 

p(x,t)y 


n- 

x 


1 


dx. 


(2.2) 
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In  Eq.  (2.2),  p(x,t)  is  the  density  the  gas,  u(x,t)  is  its  particle 
velocity,  and  p(x,t)  is  pressure.  r(j)  is  the  gamma  function.  It  assumes 
the  values  /r f  ,  1,  and  1/2  /if  for  n  =  1,2,  and  3,  respectively. 

If  one  substitutes  in  Eq.  (2.2)  for  p,  u,  and  p  the  self-similar 
solution  formulas  of  Section  2.2.  then  the  integral  is  found  to  be  inde¬ 
pendent  of  time  and  proportional  to  E0.  It  does,  however,  depend  on  n, 
y  and  K(n,y).  Therefore,  Eq.  (2.2)  provides  for  each  set  of  n  and  y  a 
corresponding  value  for  the  proportionality  factor  K(n,y).  The  integral 
has  to  be  evaluated  numerically.  We  shall  return  to  the  numerical  qua¬ 
drature  in  Section  3.2. 

The  pressure  and  the  particle  velocity  at  the  shock  are  expressed 
most  conveniently  in  terms  of  the  shock  velocity  U  =  dx  (t)/dt.  The 
latter  is,  according  to  Eq.  (2.1) 


U(t) 


_2_ 

2+n  t 


(2.3) 


The  particle  velocity  immediately  behind  the  shock  is 


us(t) 


(2.4) 


and  the  corresponding  pressure  is 


Ps(t) 


2 

Y+l 


PQ  u2(t) 


(2.5) 


where  pQ  is  the  density  of  the  ambient  gas.  The  density  immediately 
behind  the  shock  is 


Y+l 

Y-l 


(2.6) 


We  reiterate  that  Eqs.  (2.3)  through  (2.6)  are  "strong  shock"  re¬ 
lations,  i.e.,  they  are  based  on  the  assumption  that  the  ambient  pressure 
Po  is  negligible  compared  to  the  shock  pressure  p  . 

2.2.  Flow  Field  Formulas 


The  closed  form  expressions  for  the  flow  field  behind  a  strong 
blast  are  given,  e.g.,  in  References  4,  5,  6,  and  8.  In  this  report,  we 
formulate  the  solution  in  a  somewhat  simpler  but  algebraically  equivalent 
form. 
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The  formulas  express  the  distance,  particle  velocity,  pressure, 
and  density  in  terms  of  time  and  a  dimensionless  parameter  v.  The  ex¬ 
pressions  are  valid  for  ratios  of  specific  heat  that  satisfy  the  condi¬ 
tions 


1  <  y  <  7  and  y  2 


(2.7) 


The  case  y  =  2  requires  a  special  treatment  because  some  of  the  formulas 
become  singular  for  that  value  of  y.  The  solution  in  this  special  case 
is  simpler  than  the  general  solution^,  but  we  have  not  included  the 
special  solution  in  our  computer  programs. 

The  formulas  for  the  flow  variables  are 


x  =  xs(t)*y(v)  ,  (2.8) 

u  =  u  (t)  •  •  y(v)  ,  •  (2.9) 

S  V  2 

P  =  Ps(t)*g(v)  ,  (2.10) 

P  =  Ps  *  h(v)  .  (2.11) 


The  time  functions  and  p  ,  appearing  in  Eqs.  (2.8)  through  (2.11)  are 
defined  in  Section  2.1.  The  functions  y(v),  g(v),  and  h(v)  are  defined 
as  follows: 


(2.12) 


g(v) 


^ nci  (^r  fer  • 


(2.13) 


h(v)  = 


1-2C2  ^  .  C4-l 

viy-v 


C5-2C3 


7ViY  V  \  ( l~av  \ 

\v2-V! )  \viY-v2/  \l-av2y 


(2.14) 


12 


'A.  Sakurai,  " Blast  Wave  Theory" ,  Mathematics  Research  Center  Technical 
Summary  Report  497  (September  1964). 
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The  dimensionless  parameter  v  varies  between 


_  2 

Vl  (n+2) y 

► 

and 


(2.15) 


V2  (n+2)(y+l) 


The  value  v  =  Vi  corresponds  to  the  distance  x  =  0  from  the  explosion, 
and  the  value  v  =  v2  corresponds  to  a  point  on  the  shock. 


The  other  constants  in  Eqs.  (2.12)  through  (2.14)  are 


a  =  1  +  j  (Y-l)  , 


(2.16) 


Ci  = 


2+n 


C2  = 


Cq  = 


-  Y-l 


2y+n-2 


(Y-2)  2n 


Cc;  = 


y-2  5 


(2y+n-2) 2n 


(2+n) (y-l)y 


(2+n)  (ny-n+2)  (2y+n-2) (ny-n+2) 


(2+n) (ny-n+2) 


(2+n) (y-l)y 
(2-y) (ny-n+2) 


(2.17) 


The  constant  K(n,y)  enters  the  formulas  (2.8),  (2.9),  and  (2.10) 
through  the  time  functions  x  ,  u  ,  and  p  ,  respectively.  The  constant 
is  determined  by  substituting  Eql .  (2 . 8) Sthrough  (2.10)  into  the 
integral  (2.2).  After  cancelling  out  the  factor  E  one  thereby  obtains 
the  relation  ° 


1  =  K(n,y)n+2  •  B  , 
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i .  e . 


K(n,y)  =  B 


-1/ (n+2) 


(2.18) 


where 


n/2 

IT 

r£) 


<t>(v) 


(*)■■*•/ t(v)dv  ’ 

|D i+(2+n)Cl  D2/D4  +  l|  Df2  D?5+nC3  dS4  D5  , 


(2.19) 


(2.20) 


Di  =  —  , 

v2 


D2 


D*  = 


V-Vj 

V2-Vl 


1  -av 


l-av2  ’ 


D4  = 


ViY-V 
ViY-v2 

Cl 

D5  = - + 

v 


C3a 


v-Vj  1-av 


(2.21) 


Eqs .  (2.8)  through  (2.21)  provide,  together  with  the  shock  formulas 
of  Section  2.1,  explicit  expressions  for  the  flow  profiles,  i.e.,  formulas 
for  the  computation  of  the  flow  field  for  any  fixed  time  t. 

2.3.  Computation  of  Flow  Histories 

Flow  history,  i.e.,  the  description  of  the  flow  at  a  fixed  position, 
can  be  obtained  from  the  formulas  of  Sections  2.1  and  2.2  after  a  simple 
manipulation.  Let  X  be  the  distance  from  the  explosion  at  which  the  flow 
history  is  to  be  computed.  The  shock  arrival  time  at  x  =  X  is,  according 
to  Eq.  (2.1) 
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X 


(n+2) /2 


(2.22) 


t 

s 


(X) 


-(n+2)/2 

K(n,y) 


1/2 


On  the  other  hand,  one  obtains  from  Eqs.  (2.1)  and  (2.8)  the  relation 


t  =  K(n,y) 


- (n+2) /2 


(n+2) / 2 


(2.23) 


Therefore , 


t  =  ts(X)  y(v)'(n+2)/2  (2.24) 


I£  one  varies  in  Eq.  (2.24)  the  parameter  v  between  v=V2  and  v=vi, 
then  one  obtains  time  values  between  t=t  (X)  and  infinity,  i.e.,  all 
times  after  shock  arrival  at  X.  The  corresponding  values  of  the  other 
flow  variables  are  obtained  by  substituting  v  and  the  corresponding  value 
of  t  (computed  using  Eq.  (2.24))  into  the  Eqs.  (2.9),  (2.10)  and  (2.11). 

We  notice  in  passing  that  Eq.  (2.8)  can  be  simplified  for  these 
calculations  by  substituting  in  it  the  definition  of  u  (t)  and  v^.  The 
result  is 


u  =  x(t,v)  • 


(2.25) 


where  x(t,v)  is  given  by  Eq.  (2.8).  For  flow  history  computations 
x(t,v)  =  X,  of  course. 

2.4.  Particle  Path  Formulas 


Closed  formulas  for  the  particle  trajectories  (or  paths)  were  derived 
by  Laporte  and  Chang^*?.  We  shall  quote  their  formulas  in  a  somewhat 
simpler  but  algebraically  equivalent  form. 

The  trajectory  of  a  particle  which  leaves  the  shock  at  the  time 
t  =  T  is  given  by 
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t  =  T  •  w(v) 


X  -  xs(t)  •  y(v)  , 


(2.26) 


where 


„(v)  -  J- 
v2 


/YVl  -  v  \p*  /  v'vl  \p2  /  1-av  \ 
\YV1  -  v2  )  \v2-Viy  y l-av2  J 


Pj  =  (2y+n-2) (2+n) (l-y)y  •  P 


-1 


P2  =  (2+n)n(y-2) (l-y)y  •  P 


P3  =  (2+n) 2 (y-1) 2y  •  PQ  -  1  , 


(2.27) 


(2.28) 


PQ  =  2  (2ny  +  2y-n  +  2) (2-n-2y)  +  2 (ny+4)n (2-y) y  + 


(2.29) 


+2 (ny-n+2) (2+n) (y  -1) 


The  constants  v1}  v2,  and  a  are  defined  by  Eqs.  (2.15)  and  (2.16), 
respectively.  The  function  y(v)  is  defined  by  Eq.  (2.12),  and  x  (t)  is 
defined  by  Eq.  (2.1).  s 

In  Eq.  (2.26),  one  obtains  with  v=v2  the  initial  point  of  the 
particle  path,  i.e.,  a  point  on  the  shock  and  with  the  coordinates  t=T 
and  x=x  (T) .  By  decreasing  the  parameter  v  to  v=vi  one  obtains  t  and  x 
values  that  approach  infinity,  thus  covering  the  complete  particle 
trajectory. 

The  flow  variables  u,  p  and  p  are  obtained  at  any  point  of  the  path 
line  by  substituting  the  corresponding  values  of  v  and  t  into  the  Eqs. 
(2.9),  (2.10)  and  (2.11),  respectively. 

2.5.  Mach-Line  Formulas 


Closed  form  expressions  for  the  Mach-lines  in  a  strong  blast  field 
were  derived  by  Laporte  and  Chang6’7.  We  shall  present  their  formulas 
in  a  somewhat  simpler  but  algebraically  equivalent  form. 
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by 


The  Mach-lines  are  given  in  terms  of  the  dimensionless  parameter  v 


t 


ti(v)*e 


±t2(v) 


x  =  xg(t)‘  y (v)  , 


(2.30) 


where 


ti(v) 


av 

1-av 


t2(v) 


1 

— - arc  sm 

/b  j  +b2 


b i+b^ a v2 
avi (1-av) 


(2.31) 


-i  avi 

J  FT 

b2  =  [2  -  (y+1)  avj  J 


2yav 1  -  (y+1 ) 


(2.32) 


The  function  y(v)  is  defined  by  Eq.  (2.12),  x  (t)  is  defined  by  Eq.  (2.1), 
and  the  constants  Vj,  V2  and  a  are  defined  by  Eqs.  (2.15)  and  (2.16), 
respectively. 

The  constants  M+  and  M  in  Eq.  (2.30)  correspond  to  the  positive 
and  negative  sign  of  the  exponent  function  t2(v),  respectively,  and 
determine  the  two  Mach-lines  passing  through  each  point  in  the  x,t- 
plane. 

In  order  to  calculate  a  specific  characteristic  one  might  specify 
as  initial  point  either  a  point  on  the  shock  or  a  point  with  x=0  (i.e., 
at  the  location  of  the  explosion) . 

Let  the  initial  point  be  on  the  shock  at  a  distance  X  from  the 
explosion.  Then  the  corresponding  shock  arrival  time  T=t  (X)  can  be 
computed  with  Eq.  (2.22).  The  two  constants  M  and  M_  which  define  the 
two  characteristics  crossing  at  (X,T)  are  obtained  by’substituting  v=V2 
and  t=T  into  Eq.  (2.30).  The  result  is 
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T 

ti(v2) 


=  1  p+t2(v2) 

ti (v2) 


e-t2(v2) 


(2.33) 


If  the  starting  point  for  the  calculation  of  the  characteristics  is 
chosen  at  the  explosion  (x=0)  and  time  t=T,  then  the  constants  M  and  M_ 
are  obtained  from  Eq.  (2.30)  for  v=vi.  Substituting  v=v\  into  Eq.  (2.31) 
one  obtains 


t2(vl)  -  -  J  /  */,bi+b2  (2.34) 


Therefore,  in  this  case  the  formulas  for  the  constants  simplify  to 


T 

M  =  — 7 — - 
+  txCvx) 

M  =— L_ 
-  titV!) 


(2.35) 


Once  the  constants  M+  and  M_  are  found,  then  all  other  points  of 
each  characteristic  are  obtained  by  letting  the  parameter  v  vary  between 
vj  and  v2.  The  former  value  corresponds  to  a  point  with  x=0  and  the 
latter  value  corresponds  to  a  point  on  the  shock. 

The  flow  variables  u,  p,  and  p  along  the  Mach-lines  are  obtained 
by  substituting  the  corresponding  values  of  v  and  t  into  the  Eqs.  (2.9), 
(2.10),  and  (2.11),  respectively. 

3.  COMPUTER  PROGRAMS 


3.1.  General  Comments 


The  computer  program  package  consists  of  six  subroutines  for  flow 
calculations  and  two  auxiliary  routines.  One  of  the  subroutines, 

SBLPREP ,  is  a  preparation  routine.  It  computes  a  set  of  constants  from 
input  data  (ambient  conditions  and  energy  released)  and  makes  the  constants 
available  to  the  other  programs  via  a  labeled  COMMON.  The  preparation 
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routine  must  be  called  first,  before  any  of  the  other  flow  calculation 
routines  can  be  called. 

The  other  five  flow  calculation  routines  compute,  respectively,  the 
shock,  the  flow  profile  at  a  given  time,  the  flow  history  at  a  given 
distance,  a  particle  path  with  corresponding  flow  variables,  and  a  pair 
of  Mach-lines  with  corresponding  flow  variables. 

The  two  auxiliary  routines  are  used  by  the  preparation  routine  for 
its  calculation  of  the  constants. 


All  arguments  are  assumed  to  be  expressed  in  SI  base  units. 
3.2.  Preparation  Routine  SBLPREP 


The  preparation  routine  SBLPREP  computes  the  values  of  all  constants 
that  are  needed  for  flow  calculations  and  appear  in  the  formulas  of 
Section  2.  The  subroutine  must  be  called  (once)  before  any  of  the 
actual  flow  calculation  routines  can  be  used.  The  preparation  routine 
is  called  by  the  statement 

Ca 1 1  SBLPREP (N , P , TEM , GAM , AMO  L , ENCHRG , NBAD) 


The  arguments  are 


N 

P 

TEM 

GAM 

AMOL 

ENCHRG 


dimension  of  the  space  (N=l,2,  or  3); 
ambient  pressure  (in  pascals) ; 
ambient  temperature  (in  kelvins) ; 

ratio  of  specific  heats;  GAM  is  restricted  by  Eq.  (2.7); 

molar  mass  of  the  ambient  gas  (in  kg/mole) ; 

N-3 

energy  released  by  the  charge  (in  J*m  ;  one  kiloton  TNT 
equivalent  is  4.184‘10^Jj  *)  • 


NBAD  =  an  error  indicator;  it  is  set  equal  to  zero  by  SBLPREP  after 
computation  of  all  constants;  a  non- zero  NBAD  on  return 
indicates  that  the  constants  cannot  be  computed  with  the 
values  given  by  the  arguments  N  through  ENCHRG. 

^Tfie  equivalence  1  Uton  TNT  —  4. 184’ 10  J  is  given  in  reference  IS.  In 
reference  14,  p.  13,  footnote,  a  kton  TNT  is  defined  as  10 calories 
without  specifying  the  calorie  type,  and  in  the  same  reference,  pp.  13 
and  647  one  finds  the  conversion  factors  4.18,  4.187,  and  4.2’10 1%. 
According  to  reference  13,  a  calorie  has  a  value  between  4.1819J  and 
4.19002  J,  depending  on  its  type.  In  the  present  context,  the 
appropriate  type  seems  to  be  a  thermochemical  calorie"  which  is  defined 
as  4.184J. 

13  .... 

American  National  Standards  Institute,  American  So^ety  for  Tesfong 

Materials,  "Standard  for  Metric  Practice",  ASTM  E380-76,  ASTM, 

Philadelphia,  PA,  1976. 


4 Samuel  Glasstone  and  Philip  F.  Dolan,  eds,  "The  Effects  of  Nuclear 
Weapons",  3rd  edition,  U.S.  Department  of  Defense  and  U.S.  Department 
of  Energy,  1977. 
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If  any  of  the  quantities  P  through  AMOL  are  missing  (i.e.,  are  not 
positive),  then  the  following  default  values  are  used  instead  of  the 
missing  quantities-; 


P  =  101.325  kPa 

TEM  =  293.0  K 
GAM  =  1.4 

AMOL  =  28.96  g/mole 


These  values  correspond  to  a  "standard"  air.  The  use  of  default  values 
is  indicated  by  a  printed  message. 

If  the  charge  energy  ENCHRG  is  not  positive,  or  if  N  or  GAM  are  not 
within  the  indicated  ranges,  then  the  calculation  of  the  constants  is  not 
carried  out.  Instead,  a  return  with  NBAD  0  is  executed,  after  printing 
of  a  corresponding  error  message.  If  GAM  =  2.0  is  specified  by  the 
calling  routine,  then  the  computation  will  be  carried  out  with  the  default 
value  GAM  =  1.99. 

The  calculated  constants  are  stored  in  a  labeled  COMMON/SBLCOM/ 
and  made  thus  available  to  the  flow  calculation  routines.  The  contents 
of  SBLCOM  is  printed  by  the  subroutine  before  executing  a  regular  return. 

The  calculation  of  all  but  one  of  the  constants  is  trivial  and 
involves  merely  the  evaluation  of  the  lengthy  formulas  of  Section  2.  The 
computation  of  the  constant  K(n,y),  defined  by  Eqs.  (2.18)  through  (2.21), 
requires  a  numerical  quadrature,  which  is  carried  out  using  a  Romberg 
algorithm.  The  algorithm  is  implemented  by  the  auxiliary  routines 
SBLROMB  and  SBLINTE.  Because  the  integrand  <}>(v)  is  singular  at  the  lower 
limit  vi,  the  integration  is  done  in  terms  of  the  transformed  variable 


between  the  limits  zero  and  one. 

The  input  arguments  do  not  contain  the  density  of  the  ambient  air, 
which  is  needed  in  the  formulas.  The  density  is  calculated  in  SBLPREP 
from  the  given  arguments  by  the  formula. 

Pq  =  (P*AM0L) / (TEM*8 . 3143)  , 

where  8.3143  is  the  universal  gas  constant  in  J/(K*mole). 
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The  pressure  and  temperature  instead  of  density  were  chosen  as  input 
parameters  because  the  former  quantities  are  more  readily  available  when 
dealing  with  simulations  of  a  real  explosion. 

3.3.  Shock  Computation  Using  SBLSHCK 

The  subroutine  SBLSHCK  computes  the  trajectory  of  a  shock  generated 

by  a  strong  blast,  and  corresponding  shock  and  particle  velocities,  shock 

pressure,  shock  density,  and  dynamic  pressure.  The  computations  are  done 

for  a  prescribed  number  of  distances  from  the  explosion,  equidistantly 

spaced  between  X  .  and  X 

mm  max 

The  subroutine  is  called  by  the  statement 


CALL  SBLSHCK  (XMIN,  XMAX,  NRX,  X,  T,  USHCK,  P,  UPART,  RHO,  DYNPR,  NBAD) 


The  first  three  arguments  are  to  be  specified  by  the  calling  program. 
They  are 


XMIN,  XMAX  =  limits  of  distances  (in  metres)  from  the  explosion 
between  which  calculations  should  be  done. 

NRX  =  number  of  nodes  to  be  calculated  between  XMIN  and  XMAX. 

The  next  eight  arguments  (X  through  DYNPR)  are  one-dimensional  arrays 
in  which  the  calculated  results  will  be  stored: 


X  =  array  of  distances  of  the  nodes  in  metres, 

T  =  array  of  shock  arrival  times  at  the  nodes  in  seconds, 

USHCK  =  array  of  shock  velocities  (=dx/dt)  in  m/s, 

P  =  array  of  shock  pressures  in  pascals  (P  may  be  also  interpreted 
as  overpressure,  because  the  ambient  pressure  is  assumed  to  be 
negligible) ; 

UPART  =  array  of  particle  velocities  at  the  shock  in  m/s, 

RHO  =  array  of  densities  at  the  shock  in  kg/m  , 

DYNPR  =  0.5*RHO  *  UPART  **  2  =  array  of  dynamic  pressures  in  pascals. 

The  last  argument,  NBAD,  is  an  error  indicator.  It  is  set  equal  to  zero 
if  all  computations  have  been  properly  carried  out.  If  NBAD  is  not  zero 
at  return,  then  the  shock  has  not  been  calculated.  A  printed  error 
message  explains  the  reason  for  such  an  error  return.  It  can  be  caused, 
e.g.,  by  an  NRX  <0,  by  XMIN>XMAX,  etc. 
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The  formulas  for  the  calculation  of  the  shock  are  given  in  Section 


2.1. 


3.4.  Flow  Profile  Computation  Using  SBLPRQF 

The  subroutine  SBLPROF  computes  flow  profiles,  i.e.,  the  flow  field 
for  fixed  times.  The  calculations  are  based  on  formulas  for  the  flow 
field  of  a  strong  blast.  The  formulas  are  given  in  Section  2.2. 

The  subroutine  is  called  by  the  statement 


CALL  SBLPROF (T, XMIN, XMAX,NRX, X, P ,UPART, RHO,DYNPR,NBAD) 


The  subroutine  computes  for  the  time  T  a  total  of  NRX  nodes  located 
at  distances  between  XMIN  and  XMAX  from  the  explosion.  The  first  four 
arguments,  T  through  NRX,  are  set  by  the  calling  program.  They  are: 


T  =  time  (in  seconds)  after  the  explosion,  for  which  the  profile 
should  be  computed; 

XMIN,  XMAX  =  range  of  distances  (in  metres)  from  the  explosion,  between 

which  the  profile  should  be  computed  (the  actual  computation 
will  be  done  only  for  positive  distances  and  only  within 
the  blast  region); 

NRX  =  number  of  nodes  to  be  computed  between  XMIN  and  XMAX. 


The  remaining  arguments  are  set  by  the  subroutine  SBLPROF.  The 
arguments  X  through  DYNPR  are  one-dimensional  arrays  in  which  the 
subroutine  will  store  the  coordinates  of  the  calculated  nodes  and  the 
corresponding  values  of  flow  variables.  The  arguments  are: 


X  =  x-coordinates  of  the  nodes  (distances  from  the  explosion  in  metres) 
P  =  pressures  in  pascals; 

UPART  =  particle  velocities  in  m/s; 

3 

RHO  =  densities  in  kg/m  ; 

DYNPR  =  dynamic  pressures  in  pascals 
=0.5*  RHO  *  UPART  **  2 

The  last  argument,  NBAD,  is  an  error  indicator.  It  is  set  equal  to 
zero  by  the  subroutine  if  the  profile  has  been  calculated.  If  the 
profile  cannot  be  calculated  for  the  specified  arguments,  then  NBAD 
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assumes  a  non-zero  value.  In  such  a  case  also  an  error  message  is  printed 
explaining  the  reasons  for  the  error  return. 

The  computed  nodes  are  not  equidistant  in  terms  of  X.  Instead,  the 
distances  between  the  nodes  are  smaller  in  the  vicinity  of  the  shock. 

Other  arrangements  of  nodes  can  be  obtained  by  calling  the  subroutine 
repeatedly,  e.g.,  once  for  every  two  nodes,  specifying  NRX=2  and  assigning 
corresponding  values  for  XMIN  and  XMAX. 

3.5.  Flow  History  Computation  Using  SBLHIST 

The  subroutine  SBLHIST  computes  flow  histories  at  prescribed  distances 
from  the  explosion.  The  computations  are  based  on  formulas  given  in 
Section  2.3.  The  subroutine  is  called  by  the  statement 

CALL  SBLHIST (X,TMIN,TMAX,NRT,T,P, UP ART, RHO,DYNPR,NBAD) 


The  subroutine  computes  NRT  nodes  located  at  the  distance  X  from  the 
explosion  and  between  the  times  TMIN  and  TMAX  after  the  explosion.  The 
first  four  arguments,  X  through  NRT,  are  set  by  the  calling  program. 

They  are: 


X  =  distance  (in  metres)  from  the  explosion  at  which  the  history 
should  be  computed; 

TMIN, TMAX  =  minimum  and  maximum  times  (in  seconds)  after  the  explosion 
for  which  the  computations  should  be  done;  the  first  node 
will  be  calculated  for  t  =  max  (TMIN,t  ),  where  t  is  the 
shock  arrival  time  at  the  distance  X  (§q.  (2.22)); 

NRT  =  number  of  nodes  to  be  computed. 


The  remaining  arguments  are  set  by  the  subroutine.  The  arguments 
T  through  DYNPR  are  one-dimensional  arrays  in  which  the  subroutine  stores 
the  computed  nodal  values.  The  arguments  are: 


T 

P 

UPART 

RHO 

DYNPR 


array  of 
array  of 
array  of 
array  of 
array  of 


times  (in  seconds)  after  the  explosion; 
pressures  (in  pascals); 
particle  velocities  (in  m/s); 

3 

densities  (in  kg/m  ) ; 

dynamic  pressures  (in  pascals)  =  0.5*  RH0*UPART** 


2 
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The  last  argument,  NBAD,  is  an  error  indicator.  It  is  set  equal  to 
zero  by  the  subroutine  after  the  computations  have  been  completed.  In 
case  the  flow  history  cannot  be  computed  for  the  given  arguments,  NBAD 
is  assigned  a  non-zero  value  and  a  message  is  printed  explaining  the 
reason  for  the  error  return. 

The  computed  nodes  are  approximately  equidistant  in  time,  but  the 
distances  between  the  nodes  are  smaller  in  the  vicinity  of  the  shock. 

Other  node  arrangements  can  be  obtained  by  calling  the  subroutine  repeatedly 
with  proper  values  of  TMIN  and  TMAX  and  with  specifying  NRT=1  or  NRT=  2 . 

3.6.  Particle  Path  Computation  Using  SBLPATH 

The  subroutine  SBLPATH  computes  particle  trajectories  in  a  strong 
blast  flow  field.  The  initial  point  of  the  trajectory  is  a  point  on  the 
shock.  The  end  point  is  specified  by  a  maximum  time  or  a  maximum  distance. 
The  computations  are  based  on  formulas  given  in  Section  2.4. 

The  subroutine  is  called  by  the  statement 

CALL  SBLPATH ( XSHC K , TMAX , XMAX , NRNOD , X, T, P, UPART , RHO , DYNPR , NBAD ) 


The  first  four  arguments  are  set  by  the  calling  routine.  They  are: 


XSHCK  =  distance  (in  metres)  of  the  initial  point  of  the  trajectory 
from  the  explosion,  i.e.,  the  x-coordinate  of  a  shock  point; 

TMAX,  XMAX  =  coordinates  of  the  end  point  of  the  trajectory;  the  path 

will  end  either  at  time  =  TMAX,  or  at  distance  =  XMAX,  which¬ 
ever  is  reached  first; 

NRNOD  =  number  of  nodes  to  be  computed. 

The  remaining  arguments  are  set  by  the  subroutine.  The  arguments  X 
through  DYNPR  are  one-dimensional  arrays  in  which  the  subroutine  stores 
the  coordinates  and  flow  variable  values  of  the  NRNOD  computed  nodes. 

The  arguments  are 


X  = 
T  = 
P  = 
UPART  = 


RHO  = 


distances  (in  metres)  of  nodes  from  the  explosion; 

times  (in  seconds)  after  explosion; 

pressures  (in  pascals); 

particle  velocities  =  dx/dt,  (in  m/s); 

3 

densities  (in  kg/m  ) ; 
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DYNPR  =  dynamic  pressures  (in  pascals)  =  0.5*  RHO*UPART**  2. 

The  last  argument,  NBAD,  is  an  error  indicator.  It  is  set  equal  to 
zero  by  the  subroutine  if  the  particle  path  has  been  computed.  If  the 
path  cannot  be  computed  with  the  data  specified  by  the  calling  program 
then  NBAD  is  assigned  a  non- zero  value  and  a  message  is  printed  by  SBLPATH 
explaining  the  reasons  for  the  error  return. 

The  nodes  are  approximately  equidistant  in  time. 

3.7.  Mach-Line  Computation  Using  SBLMACH 

The  subroutine  SBLMACH  computes  Mach- lines  in  a  flow  field  of  a 
strong  blast.  The  user  has  to  specify  an  initial  point  of  the  Mach-lines 
either  on  the  shock,  or  at  the  site  of  the  explosion.  The  subroutine 
then  establishes  both  Mach-lines  that  pass  through  the  specified  point. 

The  actually  computed  nodes  are  located  on  segments  of  these  characteristics, 
defined  by  minimum  and  maximum  values  of  distance  and  time.  The 
computation  is  based  on  the  formulas  of  Section  2.5. 

The  subroutine  is  called  by  the  statement 
CALL  SBLMACH (XZ,TZ,XMIN,XMAX,TMIN,TMAX,NRNODE,X,T,P,UPART,RHO, DYNPR, NBAD) 


The  first  seven  arguments,  XZ  through  NRNODE ,  are  set  by  the  calling 
program.  They  are: 


XZ  =  distance  (in  metres)  of  the  initial  point  from  the  explosion 
(XZ  must  be  zero  if  TZ  is  positive  and  vice  versa) ; 

TZ  =  time  (in  seconds)  after  explosion  for  the  initial  point;  if 
TX>0  and  XZ  =  0  then  the  initial  point  is  assumed  to  be  at 
the  site  of  the  explosion;  if  TX  =  0  and  XZ  >  0  then  the 
initial  point  is  assumed  to  be  on  the  shock  at  the  distance 
XZ  from  the  explosion; 

XMIN,XMAX  _  limits  between  which  the  Mach-lines  should  be  calculated 
TMIN,TMAX  i.e.,  limits  of  a  computational  "window". 


NRNODE  =  array  containing  the  numbers  of  nodes  to  be  computed;  the  array 
contains  two  values,  NRNODE (1)  is  the  number  of  nodes  for  the 
characteristic  with  dx/dt  >  0,  and  NRNODE (2)  is  the  number  of 
nodes  for  the  other  characteristic. 


The  remaining  seven  arguments  are  set  by  the  subroutine.  The 
arguments  X  through  DYNPR  are  two-dimensional  arrays  in  which  the  sub¬ 
routine  stores  the  nodal  values  of  the  characteristics.  The  dimensions 
of  the  arrays  are  (2,NMAX),  where  NMAX  is  the  maximum  number  of  nodes 
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to  be  computed  i.e.,  NMAX  must  be  larger  than  or  equal  to  max(NRNODE(l) , 
NRN0DE(2)).  The  nodes  with  the  index  (1,  K)  belong  to  the  character¬ 
istic  with  dx/dt  >  0,  and  the  nodes  with  the  index  (2,  K)  belong  to  the 
other  characteristic. 


For  clarity,  we  show  the  arguments  with  their  dimensions.  The 
arguments  are: 


X  (2,  NMAX) 
T  (2,  NMAX) 
P  (2,  NMAX) 
UP ART (2, NMAX) 
RH0(2 ,NMAX) 
DYNPAR (2 , NMAX) 


distances  (in  metres)  of  the  nodes  from  the  explosion; 

times  (in  seconds)  after  the  explosion; 

pressures  (in  pascals) ; 

particle  velocities  (in  m/s); 

densities  (in  kg/m  ) ; 

dynamic  pressures  (in  pascals)  =  0.5*  RH0*UPART**2 


The  argument  NBAD  is  an  error  indicator.  It  is  set  equal  to  zero 
by  the  subroutine  after  the  Mach- lines  have  been  computed.  If  they  cannot 
be  computed  for  the  specified  input  arguments,  then  NBAD  is  assigned  a 
non- zero  value  and  an  error  message  is  printed  by  SBLMACH  explaining  the 
reason  for  the  error  return. 


If  only  one  Mach-line  should  be  computed,  then  the  user  may  specify 
NRNODE  =  0  for  the  other  Mach- line. 


4 .  EXAMPLE 

As  an  example  we  present  calculations  of  a  one  megaton  TNT  explosion 
in  air  assuming  spherical  symmetry.  The  calculations  were  started  by 
calling  the  preparation  subroutine  with  the  statement 

CALL  SBLPREP  (3 , 0  .  , 0 . , 0 . , 0 . , 4 . 184  El 5, NBAD) . 

The  first  argument  advises  the  program  that  the  explosion  is  spheri¬ 
cally  symmetric  ("three-dimensional") .  The  next  four  arguments  describe 
the  ambient  conditions.  Because  they  are  not  positive  the  subroutine 
SBLPREP  assigns  standard  air  default  values  for  the  ambient  pressure, 
temperature,  ratio  of  specific  heats  and  molar  mass  (see  Section  3.2). 

The  charge  energy  is  specified  as  4.184-1015J,  which  equals  one  megaton 
TNT  equivalent,  expressed  in  joules.  (The  conversion  factor  is  4.184*1(TJ 
for  one  ton  TNT  equivalent^ . ) 

The  calculation  results  are  shown  in  Figures  1  through  6.  Figures 
1  and  2  show  flow  profiles  at  0.4  seconds  after  the  explosion.  The 
profiles  were  obtained  by  calling  the  subroutine  SBLPROF .  Figures  3  and 
4  present  flow  histories  at  800  m  from  the  explosion.  They  were 
obtained  by  plotting  the  results  from  the  subroutine  SBLHIST .  Figure  5 
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PRESSURE 


Figure  1.  Pressure  Profiles  0.4  Seconds  After  the  Explosion. 

Distances  are  expressed  in  metres  and  pressures  are^expressed  in  pascals. 
Dynamic  pressure  is  defined  as  the  product  0.5  •  pu2. 
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DENSITY  VELOCITY 


\ 


Figure  2.  Velocity  and  Density  Profiles  0.4  Seconds  After  the  Explosion, 

Velocities  are  expressed  in  m/s,  densities  are  expressed  in  kg/m3,  and 
distances  are  expressed  in  metres. 
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PRESSURE 


0  2.  <4  £  S  *L0_1 


TIME 

Figure  3.  Pressure  Histories  at  800  m  from  the  Explosion. 

Times  are  expressed  in  seconds  and  pressures  are  expressed  in  pascals. 
Dynamic  pressure  is  defined  as  the  product  0.5  •  pu  . 
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VELOCITY  DENSITY 


TIME 


Figure  4.  Velocity  and  Density  Histories  at  800  m  from  the  Explosion. 

3 

Velocities  are  expressed  in  m/s,  densities  are  expressed  in  kq/m  and 
times  are  expressed  in  seconds. 
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TIME 


Figure  5.  Particle  Trajectories  of  a  one  megaton  TNT  Explosion, 
Distance  is  expressed  in  metres  and  time  is  expressed  in  seconds. 
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Figure  6.  Mach-lines  in  the  Blast  Bubble  of  a  one  megaton  TNT 
Explosion* 

Distance  is  expressed  in  metres  and  time  is  expressed  in  seconds. 
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shows  some  particle  paths  computed  using  the  subroutine  SBLPATH,  and 
Figure  6  shows  Mach-lines  computed  by  the  subroutine  SBLMACH.  The 
shock  line  in  Figure  5  and  Figure  6  was  obtained  by  calling  the 
subroutine  SBLSHCK. 

5.  LIMITS  OF  POINT  SOURCE  APPROXIMATIONS  OF  EXPLOSIONS  IN  AIR 

The  strong  blast  formulas  of  Section  2  are  based  on  assumptions 
that  are  outlined  in  Section  1  and  that  might  be  satisfied  at  distances 
not  too  close  nor  too  far  from  the  explosion.  In  order  to  make  the 
formulas  useful  as  approximations  to  real  explosions,  one  has  to  specify 
more  precisely  their  range  of  applicability.  In  this  section  we  shall 
derive  quantitative  conditons  for  which  the  strong  blast  formulas 
approximate  spherical  explosions  in  air.  The  derivations  will  be  based 
on  calculations  in  references  15  and  16.  Both  references  report  results 
obtained  by  numerical  solution  of  so-called  bursting  sphere  problems. 

A  bursting  sphere  problem  is  characterized  as  follows:  At  time 
t  =  0  one  postulates  an  ideal  gas  in  which  an  amount  of  internal  energy 
E  is  uniformly  distributed  throughout  a  sphere  of  radius  R.  The  sphere 
can  be  thought  of  as  representing  the  fireball  formed  during  the  initial 
stages  of  an  explosion  when  thermal  and  chemical  processes  dominate. 

The  time  zero  is  the  time  at  which  the  fireball  is  completely  formed  and 
the  hydrodynamic  motion  starts  by  forming  a  shock,  which  subsequently 
propagates  into  the  ambient  atmosphere.  The  initial  velocities  are 
assumed  to  be  zero  everywhere,  and  the  initial  density  is  assumed  to  be 
constant  and  equal  p,  within  the  sphere,  and  constant  and  equal  p  outside 
the  sphere.  (In  reference  15  it  is  also  assumed  that  p ^  =  pQ.)  Bet  the 
corresponding  initial  pressures  be  p1  and  pQ  ,  respectively. 

First  we  discuss  the  dependence  of  the  incident  shock  pressure  on 
distance,  as  obtained  from  strong  blast  formulas  and  flow  field 
calculations,  respectively.  The  results  from  bursting  sphere 
calculations  in  reference  15  are  displayed  in  Figures  7  and  8.  The 
curves  in  the  figures  represent  shock  pressures  in  an  ideal  gas  with 
the  ratio  of  specific  heats  y  =  1.4.  The  ordinate  in  Figure  7  is  the 
ratio  p  /p  ,  i.e.,  the  ratio  of  shock  pressure  to  initial  shock  pressure 
(at  timl  zlro) .  The  abscissa  is  the  ratio  of  distance  r  to  the  initial 
sphere  radius  R.  One  has  the  following  relations  between  the  initial 
shock  pressure  pso,  the  initial  pressure  p1  in  the  sphere  and  the  energy 


15 M.  Lutzky  and  D.  Lehto ,  "Transformations  for  Scaling  of  Close-In 
Pressures  from  Nuclear  Explosions" ,  Naval  Ordnance  Laboratory  Report 
NOLTR-66-12 ,  Silver  Spring,  MD,  March  1966. 


16Shih  Lien  Huang  and  Pei  Chi  Chou,  "Calculations  of  Expanding  Shock  Waves 
and  Late-Stage  Equivalence",  Drexel  Institute  of  Technology  Report 
125-12,  Philadelphia,  PA,  April  1968. 
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Figure  7.  Shock  Pressure  versus  Distance  for  Bursting  Spheres, 
The  curves  are  taken  from  reference  15. 
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Figure  8.  Shock  Pressure  versus  Distance  for  Bursting  Spheres 
(Sachs'  Scales), 

The  figure  is  taken  from  reference  15. 
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(5.1) 


Pso  =  0.461 *p!  =  0.0440*E/R3. 


(These  and  other  relations  in  this  section  are  derived  in  reference  15 
for  y-1.4  and  Pi/p  >>1.)  The  strength  of  the  explosion  is  characterized 
by  the  dimensionless  parameter 


The  dashed  line  in  Figure  7  represents  the  strong  blast  solution. 
According  to  Section  2.1  that  solution  is 


2  /  o  \  ^  2+n 

Ps  =  y+T  ‘  [2TK)  *  *  E  /xs 

In  the  present  case,  n  =  3,  y=1.4,  K(3,  1.4)  =  1.03278  and 

p  =  0.1567-E/r  3, 
s  s 


(5.3) 


(5.4) 


or,  by  substituting  Eq.  (5.1)  into  Eq.  (5.4), 


3.56- (R/rs)3 


(5.5) 


It  is  apparent  from  Figure  7  that  the  strong  blast  solution 
approximates  a  finite  segment  of  the  shock  curve  only  for  events  with 
S  <  0.08.  In  terms  of  initial  pressures  or  energy,  this  condition  can  be 
formulated  by  either  of  the  following  equivalent  relations 


p  /p  >'86  , 

1  so  1 o  5 

Pl/PQ  >  190  , 

E/(R3p  )  >  1950  . 
o 


(5.6) 
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If  the  event  satisfies  the  conditions  (5.6)  then  the  strong  blast 
solution  can  be  used  as  an  approximation  between  r/R  =  2.5  and  an  upper 

limit  which  depends  on  the  strength  of  the  explosion. 

The  upper  limit  for  the  approximation  can  be  obtained  from  Figure 
8  which  contains  the  same  curves  as  Figure  7,  but  uses  Sachs'  scaling 
for  its  coordinates.  Figure  8  shows  that  an  upper  limit  is  given  by 
the  condition  PS/PD  >  10* 

Both  conditions  can  be  expressed  either  in  terms  of  distances  or  in 
terms  of  pressures.  In  terms  of  distances  one  obtains  the  following 
interval  in  which  the  strong  blast  solution  can  be  used  as  an  approximation 


0.25/C 


0.25  Ve/O^pJ 


2-5  "  R  "  < 


(5.7) 


0.547  X  pi/p 


V 


V 


0.707  \  p  /p 

1  *  so  ro 


In  terms  of  pressures  one  has  the 


equivalent  conditions 


0.0101/C3 


0.0101-E/(R3Po) 
0.106-  P!/po 


>> 


ps 

Pq 


>  10. 


0.230- 


P  /P 
rso  1  o 


(5.8) 


In  order  to  check  either  of  the  conditions  one  needs  in  addition  to  the 
energy  E  also  an  estimate  of  the  fireball  radius  R.  Such  an  estimate 
is  often  not  available.  It  is,  therefore,  desirable  to  have  conditions 
which  are  independent  of  a.  If  one  is  content  with  more  restrictive 
conditions  than  (5.7)  or  (5.8),  then  one  can  replace  p  in  the  last 
Eq.  (5.8)  by  an  observed  maximum  ps  value.  Let  that  value  be  Psmax- 
Then  a  condition  for  the  applicability  of  strong  blast  shock  formulas 
to  explosions  in  air  is  the  following  set: 
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p 


s 


<  0.23 


P 


smax 


and 


(5.9) 


An  observation  of  Psmax  can  so  used  to  obtain  an  upper  limit  for 

the  fireball  radius.  That  limit  is,  according  to  Eq.  (5.1)  given  by 


R  <  0.353 


(5.10) 


The  flow  field  generated  by  a  bursting  sphere  contains  in  addition 
to  the  leading  shock  also  a  second  shock  and  a  contact  discontinuity. 
Therefore,  the  strong  blast  formulas  may  approximate  the  field  only  within 
a  strip  in  the  r,t-plane  behind  the  leading  shock,  where  the  flow  is  free 
from  these  singularities.  In  order  to  establish  limits  for  approximations 
within  the  strip  one  needs  examples  of  flows  generated  by  explosions  for  a 
reasonable  range  of  the  parameter  £.  Unfortunately  such  calculations 
have  not  been  published.  Reference  16  provides  curves  for  one  typical 
example  only.  That  example  has  the  explosion  strength  £  =  0.076,  i.e., 
it  is  a  limiting  case  where  the  strong  blast  approximation  of  the  shock 
strength  is  not  valid.  (The  calculation  is  done  for  the  density  ratio 
Pl/p  =  1.16,  instead  of  1.0,  but  this  is  probably  not  significant.) 

One  can  hypothesize  that  the  approximation  by  strong  blast  formulas  will 
be  better  for  more  powerful  explosions  than  in  this  example. 

A  comparison  between  numerical  flow  field  calculations  and  strong 
blast  formulas  is  done  in  Reference  16  in  terms  of  relative  pressure 
profiles  p/ps,  relative  particle  velocity  profiles  u/us  and  relative 
sound  speed  profiles  c/cs.  The  calculations  and  comparisons  are  presented 
for  times  corresponding  to  strong  blast  shock  positions  between  r/R  = 

3.58  and  r/R  =  7.60.  From  Figure  7  one  can  see  that  for  these  positions 
the  computed  shock  strength  differs  significantly  from  the  strong 
blast  shock  strength.  (See  curve  £  =  0.08  in  Figure  7).  Nevertheless, 
the  corresponding  relative  pressure  profiles  differ  by  no  more  than  20%, 
the  strong  blast  formulas  providing  higher  values.  (The  comparison  is 
made  within  a  strip  behind  the  leading  shock  covering  15%  of  the  blast 
bubble  radius).  The  differences  become  smaller  for  later  times.  At 
r/R  =  7.60  the  maximum  difference  is  only  about  1%. 

The  relative  particle  velocity  profiles  behave  differently.  Again, 
the  strong  blast  formulas  provide  higher  values.  The  difference  is  up 


37 


to  13%  and  the  approximation  does  not  improve  as  time  increases. 

The  relative  sound  speed  profiles  are  up  to  10%  below  the  strong 
blast  solution,  and  the  difference  does  r.ot  decrease  as  the  time  increases. 

-2 

The  gas  density  is  proportional  tc  c  and,  therefore,  the  calculated 
density  would  be  at  large  distances  up  to  20%  higher  than  predicted  by 
stro^g^last  formulas.  Consequently  the  calculated  dynamic  pressure 
u  c  )  can  be  expected  to  be  about  6%  lower  than  given  by  the 
strong  blast  formulas,  because  the  deviations  in  sound  speed  and  particle 
velocity  compensate  each  other. 

In  summary,  in  the  example,  presented  in  Reference  16,  the 
relative  incident  pressure  profiles  and  the  relative  dynamic  pressure 
profiles  can  be  reasonably  approximated  by  corresponding  strong  blast 
profiles,  particularly  at  large  distances.  This  finding  is  encouraging 
because  the  example  represents  a  case  in  which  the  initial  shock  pressure 
differs  significantly  from  the  corresponding  strong  blast  pressure.  One 
can  expect  better  approximations  for  stronger  explosions  (smaller 
parameter  c) >  but  this  hypothesis  should  be  tested  by  sample  calculations. 
Such  calculations  can  be  carried  out,  e.g.,  by  using  the  computer 
program  described  in  Reference  17. 
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LIST  OF  SYMBOLS 

a 

constant  defined  by  Eq.  (2.16) 

c 

sound  speed  (m/s) 

Cl* • ■ • ,C5 

constants  defined  by  Eq.  (2.17) 

E 

energy  contained  in  a  bursting  sphere  (J) 

Eo 

energy  released  by  the  explosion  (J,mn-3) 

g(v) 

dimensionless  function  defined  by  Eq.  (2.13) 

k(v) 

dimensionless  function  defined  by  Eq.  (2.14) 

K(n,y) 

proportionality  factor,  see  Eq.  (2.1)  and  (2.18) 

M+,  m 

factors  in  Mach  line  formulas  (2.30) 

n 

dimension  of  the  event  (n=l,2,3  for  planar,  cylindrical, 
spherical  events,  respectively) 

P 

pressure  (Pa) 

po 

ambient  pressure  (Pa) 

P1 

initial  pressure  in  a  bursting  sphere  (Pa) 

ps 

shock  pressure  (Pa) 

Pso 

initial  shock  pressure  in  a  bursting  sphere  event  (Pa) 

P 

1  smax 

maximum  shock  pressure  observed  (Pa) 

r 

distance  from  the  center  of  a  bursting  sphere 

R 

fireball  radius  (m) 

t,T 

time  after  the  explosion  (s) 

^1 (v) >  ^2  Cv)  dimensionless  functions  defined  by  Eqs.  (2.13)  and  (2.32) 


t 

s 

shock  arrival  time  (s) 

u 

particle  velocity  (m/s) 

u 

s 

particle  velocity  immediately  behind  the  shock  (m/s) 

U 

shock  velocity  (m/s) 
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LIST  OF  SYMBOLS  (Continued) 


v  running  parameter  (see  Eq.  (2.15) 

vY>v2  particular  values  of  v,  defined  by  Eq.  (2.15) 
w(v)  dimensionless  function  defined  by  Eq.  (2.27) 

x,X  distance  from  the  center  of  the  explosion  (m) 

xs,Xs  shock  distance  from  the  center  of  the  explosion  (m) 

y(v)  dimensionless  function  defined  by  Eq.  (2.12) 

y  ratio  of  specific  heats 

t,  explosion  strength  parameter  defined  by  Eq.  (5.2) 

p  gas  density  (kg/m  ) 

3 

Pq  ambient  gas  density  (kg/m  ) 

3 

p^  initial  gas  density  in  a  bursting  sphere  (kg/m  ) 

3 

Ps  gas  density  immediately  behind  the  shock  (kg/m  ) 
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APPENDIX 

LISTING  OF  COMPUTER  PROGRAMS 
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SUBROUTINE  SBLPREPUl#  AIRP#  AIRT#  AIRC# AIRM# CHARE N»N BAD) 

C  THIS  IS  PREPARATION  ROUTINE  FOR  STRONG  BLAST  COMPUTATIONS 

COHMOM/SBLCOM/NC*  GAMMA*  RO»CHEN>  V(2)  >  A*  B<2  )>  C  (  5)»  P  ( 3)  #/.K#NGOOD 

EXTERNAL  SBLINTE 

DATA(NFIRST-O) 

IF  IMFIRST.EQ.iQ)  NGOOD-O 

NFIRST-1 

NBAD-0 

IDCF-0 

PP.ES-AIRP 

if:pres.gt.o.dgoto  if 

PRES-101325.1 
IDEF-IDEF+1 
If  TEHP-AIRT 

IFITEMP.GT.O.UGOTO  25 
TEMP-293. C 
IDSF-IDEF+1 
25  GAH-AIRG 

IFI  GAM.GT. O.l) GOTO  35 

GAH-l.i4 

IDEF-IDEF+1 

35  IF [ GAM. NE. 2.0) GOTO  38 
GAM-1. 99 
IDEF-IDEF+1 
38  AMOL-AIRM 


+  *¥  *  1 

S3LPRE  2 
SBLPP.c  3 
S3LPRE  A 
S9LPRE  5 
S3LPRE  6 
SBLPRE  7 
SBLPRE  6 
SBL  PRE  9 
SBLPRcl 0 
SBL^REll 
SBLPRE12 
SBLPR61 3 
S6LPRE14 
S8LPP.E15 
SBLPP.eifc 
SBLPRE17 
SBLPRE 18 
SBLPRE19 
S3LPRE20 
SBLPRE21 
SELPRE22 
SBLPRF23 
SBLPRE24 
SBLPRE25 


IF [ AM JL.GT.O. ) GOTO  65 

AN0L-0.O2896 

IDEF-IDEF+1 

65  IF [CHAREN.CT.D. )GOTO  85 
MDAD-1 

PRINT  75.NBAD 
RETURN 

75  FORMAT ( 1H0* 10X# 30HRETURN  FROM  SBLPREP  WITH  NBAD-# 12# 
A39H  BECAUSE  CHARGE  ENERGY  IS  NOT  SPECIFIED) 

85  IF(N.GE.1.AND.N.LE.3)GOTO  101 
NBAD-2 

PRINT  95#  NBAD#N 
RETURN 


SBLPRE26 
SBLPRE27 
SBLPRt28 
SBLPRL29 
SBLPRE30 
SBLPRE31 
SBLPRE  32 
SBLPRE  33 
SBLPRE  34 
SBLPF.E35 
SBLPRE36 
SBLPRE37 
SBLPRE38 


95  FORMATdH  #  1DX*  30HRETURN  FROM  SBLPREP  WITH  NBAD-#I2# 

A1H  BECAUSE  N-#I3*24H  IS  OUTSIDE  RANGE  1  TO  3) 

101  IFtGAH,DT.H.O.AND.GAM.LT.7fcO)GOTO  105 
NBAD-3 

PRINT  103#  NBAD*  GAM 
RETURN 

103  FORMAT ( 1H0/10X* 30HRETURH  FROM  SBLPREP  WITH  NBAD-»I2* 

A  15 H  BECAUSE  GAMMA-* 1PE12.5/32H  IS  OUTSIDE  THE  RANGE  1.0  TO  7.0) 
105  NC-N 
AN-M 

GAMMA-GAM 

F.O-PP.ES*AMOL/  (8.3143*TEMP) 

CHEN-CHAREN 

V(  l)-2.0/  ( (2,I0+AII)*GAMMA) 

V ( 2) -4.0/ ( ( 2.0+AN )* (1.0+GANMA) ) 

A-1.0+AN*(GAHHA-1.0)*0.f 
C ( l)-2*0/ ( 2.0+AN) 

C ( 2 ) - ( GAMMA-1. 0 ) / ( 2. 0*G AMMA+AN-2. 0) 

C( 3)-2.D*AN*(GAMMA-2.0)/( (2.0+AM)*( AN*GAMMA-AN+2. 0) )-GAMM A* 


SBLPRE39 
SBLPRE40 
SBLPRE41 
SBLPRE42 
SBLPRE43 
SBLPRE44 
SBLPRE45 
SBLPRE46 
SBLPRE47 
SBLPRE48 
SBLPP.E49 
SBLPRE50 
SBLPRE51 
SBLPRE52 
SBLPRE  5  3 
SBLPRE54 
SBLPRE55 
SBLPRE56 
SBLPRE57 
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A  ( 2»Q+AN) * { GAHHA-1.0) / ( { 2. 0*CAMMA+aN-2. 0) *( AN*GAMMA-AM+2. 0) ) 

C { 4)*GANMA/{ GANMA-2. 0) 

C  GAHHA-2,0  IS  REPLACED  BY  1.99  AT  STATEMENT  35 

C ( 5 ) «AN* { 2.0*GAMM A+AN-2. 0 ) / { C 2.0+AN ) * ( AM*GAMMA-AN+2. 0) ) 

A+( 2. 0+AN)*GANMA*( GAMMA-1. 0)/(  2.  0*( 2. 0-GAMMA )*( AN*G AMMA-AK+2. 0 ) ) 
B(i)«(-1.0-CAMMA+2.0*CAMMA*A*V(1) )*A*V(1) /(GAMMA-1.0) 
b { 2)* {  2.l0-{  GAMMA+1.0)  *A*V  (1 ) )  /{ GAMMA-1.0) 
DA«<2.*AN)*<AN*GAMMA**2-2,0*AN*GAMMA+AN+2.0*GAMMA-2.0) 

DD«2»©*(  2. 0*AN*GAMMA+2.i0*GAMMA-AM+2.  0)*  (2.  O-AN-2.  0*GAMMA) 
A+2.iO*AM*CAMMA*<  2.0-GAMMA)  *< AN*CAMMA+4.0)  +  2.  0*  ( GAMMA+1.  0 )  *DA 
P  {  2)«=aM*GAMHA*{  2.0+AN)* ( 2.0-GAMMA)*  (GAMMA-1.0)  /D8 
P  (1)«=P  (2)  *(2.0*GAMMA+AN-2.0)/  (AN*  (GAMMA-2.0 ) ) 

P ( 3 )*P( 2) *DA/ { AN* ( 2.10-GAMMA  )*  (AM*GAMMA-AM+ 2.  0) )  - 1.  0 
C  NEXT  START  QUADRATURE  TO  COMPUTE  AK 

CALL  SE>LF.OME{  SBLINTE/  0.0#  1.0*  SF/NBaP) 

:f(nlad.£q.o)goto  135 

PP.IIIT  125/NLAD 
RETURN 

125  FORMATtlHOf  10X/  30HRETURN  FROM  SBLPP.EP  UITH  NBAD=/  14/ 

A28H  BECAUSE  OF  ERROR  IN  SBLROMB) 

135  EFtN.CQ.DFACT-i.O 

IF (N.EQ. 12)  FACT-3.:141f  9265359 
IFUI.GQ. 3)  FACT-6.28318530718 

PAK*FACT*(2«I/  (2.+AH)  )**2*(4./  (CANMA**2-1) )  *SF/  ( AN*C  (2)  ) 
AK«1.K>/DAK**(  1.0/ (2.  0+AN)  ) 

NGOOD-1 

PP.INT145/  MC/  CAM/  RO/  CHEN*  V(1  )*  V(  2)  /  A/ B  ( 1 )/ B(  2 ) /  ( C  (  J  )  /  J  =1/ 5 )/  P  ( 1)  / 
AP  ( 2 )/  P  { 3)>  AK»  NGOOD 


SBLPRE 5 

SBLPRtf 

SfeLPREi 

SLLPRct 

SBLPRE6.. 

SB l PRt63 

SBLPRE64 

Sc  LPRE65 

SBLPRE66 

SBLPRE67 

SBLFRE66 

SELPRL69 

$BLFkc7u 

SBLPRE 71 

SBLPP.E72 

SBLPRE 73 

SBLPRE74 

SELPRE75 

SBLPRB76 

SBLPRt.77 

SBLPFE76 

SBLPRE79 

SbLPRESG 

SDLPRE8 1 

SBLPRE  92 

SBLPRE 63 

SLLPRE84 

SBLPRE85 


145  FORMAT ( 1H1/ 10X/  26HC0MTENTS  OF  COMMON/SBLCOM// 

X30H  FOR  STRONG  BLAST  CALCULATIONS* //* 1H  rlOX/3HNC=/ 

AI2/  5Y.,  4HGAM=/  0PF8.5/  5X/  4HP.H0=/  1PE12. 5*  5X/  5HCHEN-/  IP  E12.  5*  /// 

Z  1H  /10X/5HV( J)«/2(2X/1PE12.5)//// 

B  1H  /10X/  2HA  =  /lPE12«i5///  1H0/10X/5HB (J)=/2(2X/1P£12«5)//// 

C  1H  * 10X/  5HC ( J) •/  5( 2X/1PE12.5 ) 

D///1H  * 10X/5  HP ( J) */3(2X/lPEI2.5)x/// 

E  IN  /10X/3HAK.=/1PE12«I5////1H  / 10X/6HNG00P=i  12/ //) 

IF(IDcF.EQ.O) GOTO  165 
PP.IMT155/  PP.ES/TEMP/GAH/  AMOL 

155  FORMAT ( 1 HO/ 10X/ 41HS0ME  OF  THE  FOLLOWING  ARE  DEFAULT  VALUES  ,/,!H 
A10X/44HASSICMED  BY  SBLPREP  AND  CORRESPONDING  TO  AIP.////1H  /10X/ 
C13HPRESSURE  =/ 1PE12. 5/ 3H  PA/// 

D1H  /10X/ 13HTE MPERATURE  -/1PE12.5/2H  K/// 

D1H  / 1GX* 13HGAMMA  =/ 1PE12.5/ // 

L1H  /1CX/13HM0LAR  MASS  =/lPE12.5/8H  KG/MOLE) 

165  RETURN 
END 


SBLPRE86 
SBLPRE87 
SBLPRE88 
SBLPRE89 
SBLPRE9G 
SBLPRc91 
SBLPRE92 
SBLPRE93 
SBLPRE94 
SELPRE95 
/ SBLPRE96 
SLLPRE97 
SDLPRE96 
SBLPPE99 
SBLPR100 
SBLPF.101 
SBLPR102 
SBLPR103 


45 


SUBR0UTII1L'  SBLROMB(F»A,BiFINT/NBAD) 

****  i 

ROMBERG  INTEGRATION  SUBROUTINE 

SDLF.OM  2 

DIMENSION  TC5/20) 

SDLP.OM  3 

HBAO-O 

SBLf.OM  4 

CALL  F(A>FA»NBAD) 

SBLROM  5 

IF(NBAO.HE.O) RETURN 

SBLP.OM  6 

CALL  F(Bj  FB»NBAD) 

SBLKOM  7 

IF(NBAD.NE.O) RETURN 

SBLF.OM  8 

T(l;l)-(FA+FD)*0.i5 

SBLROM  9 

km*  : 

SBLP0M10 

KNA-1 

SBLP.OMil 

is  o:n-floa7(kma)*z.i 

SL  LF.UU12 

F.M-0 

SLLP.0M13 

DO  2?  KA-l/KMA 

SBLR0M14 

AC* FLOAT ( 1+2* (KMA-KA) )  /DtN 

SBLROM 15 

BC*  FLOAT ( 2*KA-1 ) / DE  N 

SGLP.0M16 

ARG-AC*A+BC*B 

SBLP.0M17 

IF(B-A)  17*16*19 

SBLROMid 

17  AF:G-AMAX1(AF.G/BJ 

SBLR0M19 

ARG-ANINi(AP.G*A) 

SBLR0M20 

GOTO  20 

SCLR0M21 

18  FINT-O. 

SBLP0M22 

RETURN 

SELP.0M23 

19  ARG-ANAXKAP.G/A) 

SBLP0M24 

ARG*AMIM1(ARG*B) 

SBLR0M25 

20  CALL  F ( ARG* FM*MBAD) 

SBLF.JM26 

IF (NBAD.HE.O) RETURN 

SBLP.0M27 

FN-FN+FN 

SBLR0M28 

2S  CONTINUE 

SBLR0M29 

FM*FN/FLOAT(KMA) 

SBLR0M3G 

7(1,KN+1)*(T(1»KN)+FM)*0.5 

SBLROM 31 

THIS  IS  TRAPEZ*!  HQU  COMPUTE  ROMBERG 

SLLF.0M32 

KH-KM+1 

SLLR0M33 

KC-i 

S  BLP.0M34 

DDEM-l.i 

SBLR0M35 

3?  KC-KC+1 

SBLR0M36 

DDEN-DDEN+4.1 

SBLP.0M37 

T{KC»KM)=T<KC-l*KM)+(T(KC-l>KM)-T(KC-lf  KM-1 ) ) / < DDEM-1.  ) 

SBLK0M38 

IF (KC.L7.KN.AND.KC.lt. 15) GOTO  35 

SBLR0M39 

IF(KC.Gc.4)G0T0  45 

S6LP.0M4G 

KMA«KMA*2 

SBLR0M41 

GOTO  15 

SBLF.0M42 

45  K.AF-K.C-3 

SBLR0M43 

ITEST-0 

SELP0M44 

DO  55  KA»KAF*KC 

SBLP.0M45 

NOW  TEST  CONVERGENCE 

S  BLR0M46 

TEST=AB3(T(KA*KH)-T(KA-1*KM-1) ) 

SBLR0M47 

IF(TOST.GT.ABS(T(KC*KN) )*l.I-5  )  ITEST-1 

SBLP.0M48 

IF(TEST.LE.l. £-100) GOTO  65 

SBLR0M49 

55  CONTINUE 

S  BLR0M50 

IF  ( ITEST.iEQ.  0 )  GOTO  65 

SBLP0M51 

IF( KM.GE. 20) GOTO  65 

SBLR0M52 

KMA*KMA*2 

SBLR0M53 

GOTO  15 

SBLRUM54 

65  FIUT»T(KC*KH)*(B-A) 

SBLR0M55 

RETURN 

SBLR0M56 

END 

SBLP.0M57 
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SUBROUTINE  SBLINTE(U; F, NBAD) 

*¥**  1 

THIS  IS  THE  INTEGRAND  ROUTINE  FOR  THE  COMPUTATION  OF  AK 

SBL1NT  2 

THE  INTEGRAND  POES  HOT  CONTAIN  THE  PI-FACTOR  ETC. 

SBLINT  2 

C0HM0N/SBLC0M/N,GAM,R0,CHEN,V(2>,  A,  B(  2  >,C  <  5  ),  P  ( 3  ) ,  AK,NGOClD 

SBLINT  4 

MBAD*0 

SBLINT  5 

L'Xl*2«hC  (1)  *FLOAT  (2+N ) 

SBLINT  6 

EX2=C ( 2) *FLOAT  { N> 

SBLINT  7 

EX3*C  (5)  *2*I*C  (3}*FL0AT( M) 

SBLINT  8 

EX4-C  C4> 

SBLINT  9 

IF(VJ-l.O)  12,27,15 

S6LINT1G 

12 

I F  ( VJ )  15,25,35 

SBLINT11 

1? 

NBAD*1 

SBLINT1 2 

RETURN 

St LINT13 

25 

U=V(1) 

SBLINT14 

GOTO  37 

SBLINT15 

27 

U=V (2 ) 

SBLINT16 

GOTO  37 

SBLINT17 

35 

U*V(1)+(V(2)-V{1) )*U**(1./EX2) 

SBLINT18 

37 

D>U/V(2) 

SBLINT19 

D2=  CU-VCl) )/CVC2)— VC1>) 

3BLINT20 

D3*(l«iO-A*U)  /  (1*0-A*V(2) ) 

SBLINT2 I 

P4* (V(1>*GAN-U} /( V( 1) *GAM-V (2) ) 

SBLINT22 

P5=C (2) -( U-V {1}}*{C( 1 J/U+Ct  3) *A/ ( 1. 0-A*U) ) 

SBLINT23 

F*{ Pl**EXl*P2/D4+i«0) *D3**CX3*D4**EX4*D5 

SELINT24 

RETURN 

SBLINT25 

END 

SB LINT 2 6 
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SUBROUTINE  rCLSHCMRMIM/  F.MAX/  NR/P./T/US/P/  UP/RQ/DP/NBAD) 

C  THIS  COMPUTES  STRONG  SHOCK  BETWEEN  P.MIN  AND  RMAX 

C  HP.  «  NUMBER  OF  NODES  TO  BE  COMPUTED 

C  R  ■  DISTANCE 

C  T  «  SHOCK  ARRIVAL  TIME 

C  US  «  SHOCK  VELOCITY 

CP-  INCIDENT  SHOCK  PRESSURE 

C  UP  -  PARTICLE  VELOCITY 

C  RO  -  DENSITY 

C  DP  =  DYNAMIC  PRESSURE  -  0.5*R0*UP**2 

C  HR  AD  «  ERROR  INDICATED  BY  NBAD.NE.O 

DIMENSION  R(1)/T(1)*US(1)/P(1)/UP(1)/RQ(1)/DP(1) 

COMIIOM/SBLCOM/NC/CAM/P.OC/CHEN/ V(2 )  /  A/  B(2)  /  C(  5 )/ PC (3 ) / AK/ MG03D 

IF(NGOOD.IIE.O)  GOTO  10 

NCAD-11 

PRINT  11/NBaD 

RETURN 

21  FORMAT ( 1H0/  10X/  30Hr.ETURN  FROM  SDLSHCK  VJITH  NBAD-/I3/ 

A26H  AMD  WITHOUT  COMPUTATIONS.!/ // 1H  /10X/ 

E34HSUBR0UTIMC  SBLPP.EP  MUST  BE  CALLED  / 

C391I  BEFORE  OTHER  SBL-R0UTIN2S  CAN  BE  USED.#/) 

1C  NBAD-0 

:f(nr.ge.h)goto  25 

NBAD-1 

PRINT  15#  NBAD/NR 
RETURN 

15  FORMAT (1H0/10X/ 30HRETURN  FROM  SBLSHCK  WITH  NBAD-/I2/ 

A12H  BECAUSE  MR-/I4) 

25  R(l)-AHAXl(RHIN/0.) 

ANT-2+NC 
EX-AMT/2.  K) 

TF ACT- SORT ( ROC/CHEN) /AK**EX 

ROFACT- ( ( GAH+1.0) /( GAM-1.0) )*ROC 

DO  55  KA-l/MR 

IF ( R(KA) »,GT.O. ) GOTO  35 

T  (KA)  *0.1 

US(KA) -0 

P(KA)-0 

DP(KA)-0 

C  SINGULAR  VALUES  AT  P.-O  REPLACED  BY  ZERO 
GOTO  41 

35  T(KA)-TFACT*R(KA)**EX 
US( KA) -R(KA) / (T(KA)*EX) 

P  (KA)  -2«i0*R0C*US(  KA)**2/  ( 1,0+ GAM) 

UP (KA) -2.0*US (KA) /( 1, O+GAM) 

RO(KA)-ROFACT 

DP ( KA) «RQFACT*UP ( KA) **2*0.5 
41  IF(  KA.EQ.IIR)  GOTO  55 

P. (KA+1)  -  R(KA)  +  (  RMAX-R (1 ) )  /  FLQAT( NR-1) 

IF(R(KA+1)»GT*R(KA) )GOTO  55 
NBAD-2 

PRINT  45/MBAD/RMAX/RMIN 
RETURN 

45  FORMAT ( 1H0/ 10X/ 30HRETURN  FROM  SBLSHCK  UITH  NBAD-/I2/ 

A14H  BECAUSE  RMAX-/ 1PE12. 5/ 20H  IS  NOT  LaRGER  THAN/ 

B6H  RNIN-/1PE12.I5) 

55  CONTINUE 


i 

SBISHC  2 
SBLSHC  3 
SbLSHC  4 
SBLSHC  5 
SBLSHC  6 
SBLSHC  7 
SBLSHC  3 
SBISHC  9 
SBLSHC lw 
SLLSHC11 
SBLSHC12 
SBLSHC13 
SBLSHC14 
SBLSHC 15 
SBLSHC16 
SBLSHC 1 7 
SLLSHC16 
SBLSHC 19 
SBLSHC20 
SBLSHC21 
SBLSHC22 
SBLSHC  23 
SBLSHC24 
SBLSHC25 
SLLSHC26 
SBLSHC27 
SBLSHC  26 
SBLSHC  29 
SBLSHC  30 
SBLSHC31 
SBLSHC32 
SBLSHC  33 
SBLSHC34 
SBLSHC  3  5 
SBLSHC36 
SBLSHC  37 
SBLSHC36 
SBLSHC  39 
SBLSHC  40 
S6LSHC41 
SBLSHC42 
SELSHC43 
SBLSHC44 
SBLSHC45 
SBLSHC  46 
SBLSHC  47 
SBLSHC4  8 
SBLSHC49 
SbLSHC  50 
SBLSHC  51 
SBLSHC52 
SBLSHC  53 
SBLSHC54 
SBLSHC55 
SBLSHC56 
SBLSHC57 
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RETURN  SBl 

END  SDL 
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Cl  tft 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


subroutine  sblpruf(t#rmin#rmax#nr/R#p#up»ro#dp#nbad) 

THIS  COMPUTES  THE  FLOW  PROFILE  AT  TIME  T  BETWEEN  RMIN  AND  RMAX 

T  -  TIME  AFTER  THE  EXPLOSION 

P-MIN# P.MAX  ■  PROFILE  LIMITS 

NR  -  HUMBER  OF  MODES  TO  BE  COMPUTED 

R  -  DISTANCE  FROM  THE  EXPLOSION 

P  ■  PRESSURE 

UP  -  PARTICLE  VELOCITY 

RO  »  DENSITY 

DP  -  DYNAMIC  PRESSURE  ■  O.P*RO+UP**2 
MEAD  -  ERROR  RETURN  INDICATED  BY  NBAD.NE.O 
DIMENSION  R(l)#P(l)#UP(l)#RO(l)#DP(l) 

COilMOM/SBlCOM/NC#  GAM#  ROC,  CHEN#V(2 )»  A,  B(2)#  C  (5)#  PC<  3  )#  AK#  NGOOD 
Y(D)»(V(2)/D)**C(1)*((D-V(1)) /(V(2)-V(l)) )**C(2) 

A  *(  (l.!-A*D)/(l.l-A*V(2))  )**C(3) 

G ( D )■ ( D/V(2) ) **(C(l)*FLOAT(MC) )*( (V(l)*GAM-D)/( V ( 1 ) *GAM-V ( 2 ) ) 
A  4)*((l,hA*D)/(l,hA*V(2)) )**(2,*C(5) J 
H(D) * ( ( D-V(l)  )/CV(2)-V(l)  ))**(1.«-2..*C(2)) 

A  *( (V(1)*GAM-D)/(V(1)*GAM-V(2) ) )** (C ( 4)-l. ) 

B  *( (l.^A*D)/(l.hA*V(2)) )**(2.*C(5)-2.*C(3) ) 
U(D)*VC1)+(V(2)-V(1) )*D**(1./C(2) ) 
if(ngood.me.o)  goto  10 

MBAD-11 
PRINT  11#  MB  AD 
RETURN 

II  FORMAT ( 1HO#  10X#  30HRETUP.M  FROM  SBLPP.OF  WITH  NBA0-»I3# 

A26H  AND  WITHOUT  COMPUTATIONS.!# /#  1H  #10X# 

B34HSUBR0UTINE  SBLPREP  MUST  BE  CALLED  # 

C39H  BEFORE  OTHER  SBL-ROUTINES  CAN  BE  USED.#/) 

10  NBAD-0 

IF(T.GT.O,[)GOTO  25 
NBAD-1 

PRINT  15»NBAD#T 
RETURN 

15  FORMAT ( 1H0#  10X#  30HP.ETURN  FROM  SBLPP.OF  WITH  NBAD-#I2# 

A11H  BECAUSE  T-#1PE12.5) 

25  IF(NR.GE.a)GOTO  45 
NBAD-2 

PRINT  35#NBAD#NR 
RETURN 

35  FORMATt  1H0#  10X#  30HP.ETUF.N  FROM  SBl.PROF  WITH  NBAD«#I2# 

A12H  BECAUSE  NR»#I4) 

45  AHT-MC+2 

RS»  AK*(C'IEN/R  OC  )**(!. /ANT  )*T**(  2. /ANT) 

USHCK-2.0*P.S/(T*ANT) 

US»2.0*USHCK/ (l.+GAM) 

PS*2.D*R0C*USHCK**2/  ( l.'+GAM ) 
r.OS-( GAN+l.i)  *ROC/ (GAM-1.  ) 

F;NI«AMAX1(  RMIN#  0.1) 

RMI*AMIM1(P.MI#RS) 

P.MA-AMIM1  (RMAX# RS) 

NEXT  FIND  PARAMETER  V  CORRESPONDING  TO  RMI 
Xl-O.C 
Fl-O.i 
X2-1.0 
F2-1.I0 
ICNT-0 


+  ¥**  1 
SBLPRO  2 
SBLPRO  3 
SbLPRO  4 
SBLPRO  5 
SBLPRO  6 
SBLPRO  7 
SBLPRO  6 
SBLPRO  9 
SBLPR01C 
SBLPR011 
SBLPR012 
SBLPRQ13 
SBLPR014 
SBLPR015 
) **C ( SBLPR016 
SBLPR017 
SBLPR018 
SBLPR019 
SBLPR02G 
SBLPRO?.! 
SCLPR022 
SBLPR023 
SBLPR024 
SBLPR025 
SELPR026 
SBLPR027 
SBLPR028 
SBLPR029 
SBLPR030 
SBLPR031 
SBLPR032 
SBLPR033 
SCLPR034 
SBLPR035 
SBLPR036 
SBLPP.037 
SBLPR038 
SBLPR039 
SBLPR040 
SBLPR041 
SBLPR042 
SBLPR043 
SBLPR044 
SBLPR045 
SBLPR046 
SBLPR047 
SBLPR048 
SBLPR049 
SBLPR050 
SBLPP.O  51 
SBLPR052 
SBLPR053 
SBLPR054 
SBLPP.055 
SBLPR056 
SBLPR057 
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:f{rmi.gt.o.')coto  55 

SBLPRG58 

vni«o*o 

SDL  PRO 5 9 

GOTO  85 

SELPR060 

55 

:F(RMI«1T»P.S)G0T0  65 

SB  L  PP.061 

VNI-1.0 

SoLPR0b2 

GOTO  35 

S6LPRU63 

65 

>:3«X2-(X2-X1)*(F2-RMI/RS)/<F2-F1) 

SBLPR064 

X’J«U(X3) 

SGLPR065 

F3«Y(XVJ) 

SBl PR066 

IF  { ABS  { F3-RMI  /RS)*LT#1« 0E-7*F.MI/RS )  GOTO  75 

SBLPRU67 

ICNT-ICMT+l 

SBLPP.066 

:fccnt.gt»£0)goto  75 

SBLPRU69 

:F{F3-P.MI/RSJ  56*75,  56 

SBLPR070 

56 

>!>X3 

S8LPR071 

FI*  F3 

SBLPR072 

GOTO  65 

SBlPP.073 

58 

N2-X3 

SBLPR074 

F2*F3 

SBLPR075 

GOTO  65 

SBLPF076 

75 

VHI*X3 

SE.LPR077 

85 

R(1)«F;S*V(U(VH:>) 

SbLPR076 

F’(1)«PS*G(U(VMI)) 

SBLPRU79 

UP ( 1)  =US* (R{  1 ) /RS  )*(  VMV'MI )  /V  (2) ) 

SBLPP08C 

RO(l)-RUS*H(U(VMI)> 

SbLPR081 

DP ( 1) *RO{ 1)*UP ( 1) **2*0«  5 

SBLPR062 

IF  (MR*EQ*.l)  RETURN 

SBL PRO 83 

IF (RMA«GT*RMI )GOTO  105 

SLLPP.084 

I1DhD*3 

SBLPR085 

PRINT  95/ NBAD/RMAX 

SBLPR086 

RETURN 

SBLPR067 

95 

F ORMAT { 1H0/  10X,  30HP.ETUP.M  FROM  SBLPP.OF 

WITH  NBAD-,12, 

SBLPR086 

h 14H  BECAUSD  RMAX=, 1PE12. 5, 18H  IS  OUTSIDE  RANGE ) 

SBLPR069 

135 

I F ( RMA»  LT«RS ) GOTO  115 

S  BLPR090 

VHA=1«00 

SBLPP.091 

GOTO  145 

SBLPR092 

NOU  COMPUTE  PARAMETER  VMA  CORRESPONDING 

TO  P.MA 

SBLPR093 

115 

ICMT-0 

SBLPR094 

Xl-VMI 

SBLPP.095 

Fl*Y(lJ{Xi) ) 

SfcLPR096 

X2-1..0 

SBLPR097 

F2=1»I0 

S8LPR096 

125 

N3=X2-{  X2-X1 )  *{  F2-RMA/P.S )  /  (  F2-F1) 

SE.LPR099 

XIJ-W(X3) 

SBLPP.100 

F3=Y(XW) 

SbLPRIOl 

:F(ADS{F3-RMA/RS).LT.1.0E-7  *RMA/RS) 

GOTO  135 

SbLPP.102 

ICMT=:CNT+1 

SLLPR1G3 

:f(icnt.gt.20)goto  135 

SBLPR104 

IF(F3-RMA/RS)126/135/120 

SBLPR105 

126 

X1*X3 

SbLPklG6 

F1  =  F3 

SBLPR107 

GOTO  125 

SBLPR1C8 

126 

X  2*X3 

SBLPR109 

F2=F3 

SBLPR11G 

GOTO  125 

SBLPR111 

135 

VMA-X3 

SBLFP112 

145 

DO  155  KA«2,  NR 

SBLPR113 

VK= VMI+ { VMA-VMI  )*FLOAT(  KA-1 )/  FLOAT  ( NP. 

-1) 

SEtLPP.114 
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VK-U(VK) 

F.(KA)«RS*Y(l'K) 

P(KA) =PC*G(VK) 

UP(KA)-US*(r;(KA)/RS)*(VK/V(2)  ) 
P.O(KA)  =RO£*H  ( VK ) 

PP  ( KA )  -P.0 { KA )  *UP ( KA )  **2*0. 5 
155  COilTIIIUE 

return 

cud 


SBLPR115 

SBLPR116 

SbLPK117 

SBLPR118 

SBLPR119 

S&IPR120 

SBLPR121 

SBLPR122 

SBLPR123 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE 
THIS  COMPUTES 


R  ■  DISTANCE 


■  HISTORY 
OF  NODES 


LIMITS 

TO  BE  COMPUTED 


CBLHIST(P.*  THIN# T MAX#  NR#  T#  P#  UP*  RO#  DP#  NBAD ) 

THE  FLOW  HISTORY  AT  DISTANCE  P.  AND  TIME  BETWtEN  TMIN 

wTAN''" 

TMIfl#  TMAX 
HP.  -  NUMBER 
T  -  TIME 
P  -  PRESSURE 
'JP  -  PARTICLE  VELOCITY 
RO  -  DENSITY 

DP  -  DYNAMIC  PRESSURE  «  0.5*R0*UP**2 
NBAD  -  ERROR  RETURN  IS  INDICATED  BY  NBAD.NC.O 
DIMENSION  T(1)*P( l)*UP(l)*RO( 1)#DP( 1) 

COMMON/ SB LC ON /NC* CAM* ROC* CHEN# V( 2 )#  A*B(2)*C(5)»PC(3)# AK#NGCQD 
Y(D)-(V(2)/D)**C(1)*((D-V(1))/(V(2)-V(D)  )**C(2) 

A  *(  (l.i-D*A)/(l.HV(2)*A)  )**C(3) 

C ( D ) ■ ( D / V ( 2 ) ) ** ( C ( 1 ) *  F  L  OAT ( NC } ) * ( ( V ( 1 ) +GAM-D ) /  ( V ( 1 ) *G AM-V ( 2 ) ) ) * 
A  4)*((l«-D*A}/(l.h-V(2)*A)  )**<  2.  *C  (5)) 

H(D) - ( ( D-V(l)  )  /  ( V  (2  )-V(  1)  ))**{  1.-2.  *C (2 ) ) 

A  *( ( V( 1) *GAM-D) /( V( 1) *GAM-V (2) ) ) ** ( C ( A) -1. ) 

B  *(  (l.HD*A)  /(1»M/(2)*A)  ) **(  2«  *C  (5  )-2.  *C  (3) ) 

U ( D ) ■ V ( 1 } + ( V ( 2 } -V ( 1 )  )*D**(1./C(2)  ) 

IF(NGDOD.NE.O)  GOTO  10 
MBAD-11 
PRINT  11*  IIBAD 
RETURN 

11  FORMAT ( 1HG*  10X*  30HRETUF.N  FROM  SBLHIST  UITH  NBAD 
A26H  AMD  WITHOUT  COMPUTATIONS.!/ // 1H  *10X* 

B34HSUBP.0UTIME  SBLPREP  MUST  BL  CALLED  * 

C39H  BEFORE  OTHER  SBL-ROUTINES  CAN  BE  USED.*/) 

10  NBAD-0 

IF(  R.GT.lO.I)  COTO  25 
MBAD-1 

PRINT  15# NBAD* R 
RETURN 

rDRNATdHOi  10X*  30HRETURN  FP0I1  SBLHIST  UITH  NBAD*=»I2* 

n  r 


*TMS3LHIS 


Li* 


15 


25 


A11H  BECAUSE  R*#1PC12.5) 
27 


:f(nr.ce.ii)goto 

NBAD-2 

PRINT  26* NBAD* NR 
RETURN 

26  FORMAT ( 1H0*  10X*  30HF.ETURN  FROM  SBLHIST  WITH  NBAD-,  12* 
A12H  BECAUSE  NR-*I4) 

27  ANT-0.5*FL0AT(2+NC) 

ROS-( GAN+1. ) *ROC/ (GAM-1.  ) 

TS-SQRT (ROC/ CHEN) *(R/AK)**ANT 

THIS  IS  SHOCK  ARRIVAL  TIME  AT  R 
TMI-ANAX1 (TMIN* TS ) 

IF ( TNI • GT»  TS) GOTO  35 

VMI-i.O 

COTO  65 

NEXT  COMPUTE  PARAMETER  VMI  CORRESPONDING  TO  TMI 
35  ICNT-0 

PF-(T3/THI)**(i./AH7) 

Xl-0.0 

Fl-0 

Y.  2  -1.0 

F2-1.I 


SBLHIS 
S6LHIS 
SBLHIS 
S3LHIS 
SBLHIS 
SbLHiS 
SBLHIS 
SBLHIS1C 
SCLHIS11 
SBLUIS12 
SBLHIS13 
SBlHISIA 
SBLHIS15 
*0  (SBLHIS16 
SBLHIS17 
SBLH1S18 
SBLHIS19 
S6LHIS20 
SBLHIS21 
SBLHIS22 
SBLH1S23 
SBLHIS24 
SBLHIS25 
SBLHIS26 
SBLHIS27 
SBLHIS2B 
SELHIS29 
SBLHIS30 
SBLHIS31 
SBLHIS32 
SBLHIS33 
SBLHIS34 
SBLHIS35 
SBLHIS36 
SBLHIS37 
SBLHIS36 
S6LHIS39 
S6LHIS40 
SBLHIS41 
SBLHIS42 
SBLHIS43 
SBL HIS44 
S6LHIS45 
SBLHIS46 
SBLHIS47 
SBLHIS46 
SBLHIS49 
S6LHIS5C 
SBLHIS51 
SBLHIS52 
SCLHI253 
SBLHIS54 
SBLHIS55 
SBLHIS56 
SBLHIS57 
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X3*X2-{ X2-X1) *{ F2-DF) /(F2-F1) 

XU-H(X3) 

F3*Y( XW) 

IF(ABS( F3-DF) •LT«l«0E-7  *DF)GOTO  55 
ICNT-ICNT+1 
IF(ICNT.GT.£0)G0T0  55 
IF ( F3-DF ) 46/55/ 48 
46  X1-X3 
F1-F3 
GOTO  45 
4C  X2-X3 
F2-F3 
GOTO  45 
55  VMI-X3 
65  XVJ-W(VHI) 

YKA-Y(XM) 

T(1)«TS*YKA**(-ANT> 

USHCK* (R/YKA) / (T( 1) *ANT) 

P  { 1) *2#i*{  ROC  /  (l.+GAM)  )*USHCK**2+G (XW) 

UP(1)-R*XU/T<1) 

RO(  1)  *P.OS*H(  XW) 

DP(l)«RO(l)*UP(l) **2*0.5 
IF (NR.EQ.l) RETURN 
IF(TMAX.GT.TNI)GOTO  85 
I1BAD-3 

PRINT  75/ NBAD/TMAX 
RETURN 

75  FORMAT ( 1H0, 10X, 30HRETURN  FROM  SGLHIST  WITH  NBAD«, 14, 
A14H  BECAUSE  TMAX-, 1P212.5/17H  IS  OUTSIDE  RANGE) 

35  ICNT-G 

DF" (TS/THAX) ** ( 1* /ANT ) 

C  NOW  FIND  VMA  CORRESPONDING  TO  TIME«TMAX 
Xl-VHI 

F1»Y(U(VHI) ) 

X2-0.O 

F2«0.I0 

65  X3*X2-{X2“X1 ) *{ F2-DF ) / (F2-F1) 

XU»W(X3) 

F3*Y(XW) 

I F  C  AB S  (  F3-DF)  •  LT»iDF*l«0E-7  )GOTO  105 
ICNT-ICNT+l 

IF(ICNT.GT.20)G0T0  105 
IF  ( F3-DF )  96»  105,98 
96  X2*X3 
F2*F3 
GOTO  95 
98  X1«X3 
F1-F3 
GOTO  95 
105  VMA-X3 

DO  115  KA«2,  NR 

VK.*VMI+(  VMA-VMT  )*FL0AT(KA-1)/FL0AT(  NR-1) 

XW-U(VK) 

YKA«Y(XU) 

T(KA)«TS*YKA**(-ANT) 

RS-R/YKA 

USHCK«RS/(T(KA)*ANT) 


SBLHIS58 

SBLHIS59 

SBLHIS60 

SBLHIS61 

SBLHIS62 

SBLHIS63 

SBLIIIS64 

SBLHIS65 

SBLHIS6& 

SBLHIS67 

SBLHIS6e 

SBLHIS69 

SBLHIS7C 

SBLHIS71 

SBLHIS72 

SBLHIS73 

SELHIS74 

S6LUIS75 

SBLHIS76 

SBLHIS77 

SBLHIS78 

SBLHIS79 

SBLHIS80 

SBLHIS81 

SBLI1IS82 

SBLHIS83 

SBLHIS84 

SBLHIS85 

SBLHIS86 

SBLHIS87 

SBLHIS88 

SBLHIS89 

SBLHIS90 

SBLHIS91 

SBLH1S92 

SBLHIS93 

SBLHIS94 

SBLHIS95 

SBLHIS96 

SBLHIS97 

SBLHIS98 

SBLHIS99 

SELHI100 

SBLHI101 

SBLHI102 

SBLHI103 

SBLHI1C4 

SELHI105 

SBLHI106 

SBLHI107 

SBLHI108 

SBLHI109 

SBLHI110 

SBLHI111 

S6LHI112 

SBLHI113 

SBLHI114 
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P(KA)=2.*(R0Cm.:+GAM))*USHCK**2*G(XW)  S6LHI115 

UP  (  KA) *R*XU/T  (KA)  ^BLHIllb 

RO  ( KA) =ROS*H  ( XVI )  SBLHI117 

DP  ( Y.t. )  «RO  ( KA )  *UP  ( KA )  **2*0 .  f>  SBL  HI  1 1 8 

115  COMTIMUc  SBLHI119 

R2TURH  SBLHI120 

END  SELHI121 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE  SDLPATHCRZ/P.NAX/ TMAX» NR#  Rz  T/PiUP/RUzDP/NbAC  )  **** 

THIS  COMPUTES  A  PARTICLE  PATH  STARTING  ON  THE  SHOCK  AT  k-RZ  SBLPAT 

RZ  *  DISTAMCc  OF  FIRST  NODE  FROM  EXPLOSION/  LOCATED  ON  THE  SHOCK  SBLPAT 

RI1AX/TMAX  «  PATH  ENDS  UHEN  EITHER  OF  THESE  IS  REACHED  SBLPAT 

IIP.  «  DUMBER  OF  NODES  TO  BE  COMPUTED  SBLPAT 

R  *  DISTAMCE  FROM  THE  EXPLOSION  SBLPAT 

T  -  TIME  AFTER  EXPLOSION  SBLPAT 

P  -  PRESSURE  SBLPAT 

UP  *  PARTICLE  VELOCITY  SBLPAT 

P.O  «  DEMSITY  SBLPAT1G 

DP  «  DYNAMIC  PRESSURE  -  0.f-*R0*UP**2  SBLPATil 

MBAD  «  ERROR  RETURN  WILL  BE  INDICATED  BY  f’EAD.NE.O  SBLPAT12 

DIMENSION  R(l)/T(l)/P(l)/UP(L)zRO(l)/DP(l)  SBLPAT13 

COMMJN/SBLCON/MCz GAM/  ROC*  CHEN*  V(2)z  Az B(2)z C <5)z PC (3)z AK> NGOOD  S6LPAT14 
Y (D)- <  V( 2 ) /P) **C ( 1) * ( (D-V (1))/(V(2)-V(1)))**C(2)  SBLPAT15 

A  *(  ( 1*“D*A )  /  ( l*h-VC  2  )*A)  )**C  (3 )  SBLPAT16 

G(D) = (D/V( 2) ) **(C (1 )* FLOAT CMC ) )*( ( V ( 1 )*CAM-D ) / ( V ( 1} *GaM-V (2) ) )**C( SBLPATi? 


A  4)*(  ( l.-D*A)/(  1«K'(2)*A)  )**(2.*C  (5)) 

H(D)*  ( ( D-VC 1)  ) / (V (2)-V( 1)  ) )  ** (l«-2«  *C (2) ) 

A  *<  <V(1)*GAM-D)/<V<1)*GAM-V(2))  )** ( C < 4) -1. ) 

B  *(  C1.!-D*A)  / (i«K’(2 )*A)  )**(2.*C(5)-2»*C(3)  ) 

UCD) ■ (D/VC2) ) *( (V(l)*GAM-D)/( V(1)*GAM-Y(2) ) )**PC<1) 

A  *(  (D-V(l) )  /  (V(2)-V(l) )  )**PC( 2)*(  (i;-D*A)  /  ( l.^-VC  2)*A)  )**PC  (  3) 
P.SF  (D)«AK*( CHEN/ROC )**< l.;/FL0AT(2+NC )  )*D«*(  2, /FL0ATC2+NC)  ) 
X(D)-V(1)+(V(2)-V(1))  *D**(-l«iyPC  (  2) ) 

IF( NGOOD. NE.O)  GOTO  10 
NDAD-11 
PRINT  11/ MBAD 
RETURN 

11  FORMAT ( 1H0/  10X/  30HRETURN  FROM  SBLPATH  IJITH  NBA0-,I3/ 

A26H  AMD  WITHOUT  COMPUTATIONS.!. /i  1H  /10X/ 


SBLPAT16 

SBLPAT19 

SBLPAT20 

SBLPAT21 

SBLPAT22 

SBLPAT23 

S6LPAT24 

SBLPaT25 

SBLPAT26 

SBLPAT27 

SBLPAT26 

SBLPAT29 

SBLPAT3C 

SBLPAT31 


B34HSUBR0UTINE  SBLPREP  MUST  BE  CALLED  / 

C39H  BEFORE  OTHER  SBL-P.OUTIMES  CAN  BE  USED.//) 

10  NCAD-0 

IF ( RZ.GT.O.l)  GOTO  2f 
MBAD-1 

PRINT  15/NBAD/RZ 
RETURN 

If  TORHATCIHO/IOX/ 30HP.ETURN  FROM  SBLPATH  WITH  NBAD-z  12/ 
A12H  BECAUSE  RZ-/1PL12.5) 

25  IFCNR.GT.OGOTO  45 
NDAD-2 

PRINT  35/NBAD/NR 
RETURN 

35  F0RI1ATC  lHOz  10X/  30HRETURN  FROM  SBLPATH  WITH  NBAD-/ 12/ 
A12H  BECAUSE  NR-/ 14) 

45  RCD-RZ 
ANT-2+NC 

T  (1)=AK**  (-0.5*  ANT)  *SQP.T(  ROC/CHEN  )*RZ**(0.f  *ANT) 
USHCK-C  2.10/ANT)  *RZ/T  ( 1) 

P(l)  =  (2.O/(GAM+l»l))*R0C*USHCK**2 
UP ( 1) =( 2.0/C  GAH+1. ) ) +USHCK 
ROC  1)  *  ( ( GAM+1,  )  /  ( GAH-1,  ) )  *P,OC 
DP(1)=R0(1)*UP(1)**2*0.5 
IF (NR.EQ.l) RETURN 

IF  ( TMAX.6T..T  ( 1 ) .  AND.  RMAX.  GT.  P.  ( 1 ) )  GOTO  65 
NBAD-3 


SBLPAT32 
5BLPAT33 
SBLPAT34 
SBLPAT35 
SBLPAT36 
SBLPAT37 
SBLPAT36 
SBLPAT39 
SBLPAT  40 
SBLPAT41 
SBLPAT42 
SBLPAT43 
S6LPAT44 
SBl PAT45 
SBLPAT  46 
SBLPAT  47 
SBLPAT48 
SBLPAT49 
SBLPAT50 
SBLPAT51 
SELPAT52 
SBLPAT53 
SBLPAT54 
SBLPAT55 
SBLPAT56 
SBLPAT57 
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PRINT  55»  NBAD*  TMAXi  T  ( 1)  j  RMAX,  Rt  1) 

RETURN 

55  FORMAT ( 1MO/  10X*  30HRETURN  FROM  SBLFATH  WITH  NBAD=iI2j> 

A14H  BECAUSE  TMAX«.rlPE12.5,13H  .  LE.  TSHOCK  ■*»  1P£12.5» /» 

D1H  /10X/16H0R  BECAUSE  RMaX=.»1PE12.5>9H  .Lt«  RZ«,1PE12.5) 
65  XCHT-9 

C  MOW  COMPUTE  PARAMETER  YHA  CORRESPONDING  TO  TMAX 
DT-TCD/TMAX 
Xl«O.IO 
F1=0*0 
X2-1.0 
F  2*  1.  iO 

75  X3=X2-( X2-X1 ) *( F2-DT ) /( F2-F1) 

F3«1.!/U(X(X3)  ) 

IF(ABS(F3-DT)  .LT.Il.iOE-7*DT)  GOTO  65 

ICMT-ICMT+1 

IF ( ICNT.GT. 20JGOTO  85 

IF(F3-DT)  76»  85>78 

76  X1=X3 
F>F3 
GOTO  75 

78  M2=X3 
F2=F3 
GOTO  75 
55  VTMA»X3 

RTNAX*RSF { TMAX) *Y (X( VTHAX ) ) 

IF( RTMAX.LT.  RMAX) GOTO  105 
C  BRANCH  IF  EMD  OF  PATH  DETERMINED  BY  TMAX 
C  ELSE  COMPUTE  PARAMETER  CORRESPOMDIMG  TO  RMAX 
ICMT=0 
Xi-VTHA): 

Fl-RTMAX 

X2=1.0 

F2-R2 

95  X3=X2-( X2-X1) *( F 2-RMAX) / <  F2-F 1) 

DT=T( 1 ) *W (X ( X3) ) 

F3=r.SF(D7)*Y(X<X3)) 

:F(ABS(F3-RMAX).LT.il.0E-7*RMAX)  GOTO  105 
ICNT-ICNT+1 

IF  ( ICMT.iGTt20)GOT0  105 
IF( F3-RMAX)  96> 105»  98 

96  X2=X3 
F  2=F3 
C0T0  95 

98  X1=X3 
F1=F3 
GOTO  95 

105  YMAX-1./X3 
YMIU-1.0 
DO  135  KA  =  2»  MR 

YK= YMIM+ ( YMAX-YMIM) AFLOAT (KA-1) /FLOAT (  MR-1) 

VK=V( 1) +( V( 2) -V (1 ) ) *YK**( 1. /PC ( 2) ) 

T  (KA)  =T(l)*l.'(VK) 

RS=RSF(T( KA) ) 

R(KA)=RS*Y(VK) 

USHCK=( RS/T ( KA) )*2.©/FL0AT( 2+NC ) 

P (KA) * ( 2. 0/ ( l.+GAM) ) *R0C*USHCK*+2*G (VK ) 


S3LPAT58 
S3LPAT59 
SBLPAT6C 
SBLPAT61 
SBLPAT62 
SBLPAT63 
SBLPAT64 
SELPAT  65 
SBLPAT66 
SBLPA767 
SBLPAT68 
SSI  PAT 69 
SBt  PAT70 
SSLPAT71 
SBLPAT72 
SBLPAT73 
SSLPAT74 
SBLPAT75 
SBLPAT76 
SBLPAT77 
SBLPAT78 
SBLPAT79 
SBLPAT80 
SBLPAT81 
SBLPAT82 
SCLPAT83 
SBLPAT84 
SBl PAT85 
SBLRAT86 
SBLPAT87 
SBLPAT88 
SBLPAT89 
SBLPAT90 
SBLPAT91 
SBLPAT92 
SBLPAT93 
SBLPAT94 
SBLPAT95 
SBLPAT96 
SBLPAT97 
SBLPAT98 
SBLPAT99 
S  BLPA100 
SBLPA-C1 
SBLPA1G2 
SBLPA103 
SBLPA104 
SBLPA105 
SBLPA106 
SBLPA1C7 
SBLPA108 
SBLPA109 
SBLPAL10 
SBLPA111 
SBLPA112 
SBLPA113 
SBLPA114 
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UPCKA)"RCKA)*VK/T  CK.A)  SBLPA115 

P.0(  KA)  ■  ( ( GAM+l.i) /  C GAM-1* )  )*P.OC*H(  VK )  SBLPA116 

DP (KA)-RO(KA)*UP<KA) **2*0.5  SBLPAU7 

135  COMTIMUE  ‘‘BLPA118 

RETURN  SBtPAil1* 

END  SBLPA120 
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o  o  o  o  o  o  r>  o  o  o  r>  o  o  o  o  o 


c 


SUBROUTINE  SDLMACH(RZ*  TZ*  RMIllj  RMAX*  TMIN*  TMAX*  NR*  ****  1 

A  R*T*P*UP*R0*DP.'I1EAU)  SbLMAC  ? 


THIS  COMPUTES  MACH-LINES  AMD  CORRESPONDING  FLOW  FIELD 
RZ*TZ  «  INITIAL  POINT  OF  THE  MACH  LINES. i  EITHER  RZ  OR  TZ  MUST 
BE  ZERO. I  IF  RZ-0*  THEN  COMPUTATION  HILL  START 
AT  THE  EXPLOSION  AND  TINE  TZ. 

IF  TZ-0#  THEN  COMPUTATION  UILL  START  ON  THt 
SHOCK  AT  DISTANCE  RZ  AND  CORRESPONDING  TIME 

RMIN* - * TMAX  =  LIMITS  FOR  THE  COMPUTATIONS 

MR ( 2 )  =  MUMPER  OF  NODES  TO  BE  COMPUTED. 

NR ( 1 )  CORRESPONDS  TO  CHARACTERISTIC  WITH  DT/DR.GT.O 
R(2j i)  =  DISTANCE  FROM  THE  EXPLOSION 
7(2, 1)  =  TIME  FROM  THE  EXPLOSION 

P  ( 2*  1 )  =  PRESSURE 

UP(2*i)  -  PARTICLE  VELOCITY 
R0<2*1)  =  DENSITY 

DP (2* 1 )  =  DYNAMIC  PRESSURE  =  U.5*P.O*UP**2 

MB AD  -  ERROR  RETURN  WILL  BE  INDICATED  BY  NBAD.NE.O 

DIMENSION  HR(2)*R<2*1)*T<2*1)*P<2*1)*R0(2*1)*DP(2*1)*UP<2*1) 
COMMOM/SBLCOM/NC*  GAM*  ROC* CHEN* Y( 2)*  A*  B( 2) *C ( 5) »  PC ( 3 )*  AK* NGQOD 
Y(D)-(V(2)/D)**C(i)*((D-V(l>>/(V(2)-V(l)) )**C(2) 

A  *(  (l.l-D*A)  /{ l.h-V(  2)*A)  )**C  (3) 

G(D)=(D/V(2))**(C(1)*FL0AT(NC ) )* ( ( V (1 )*GAM-D ) / ( V ( 1 ) *GAM-V (2 ) ) 
A  4)*((i.-D*A)/(l,t-V(2)*A))**<2.*C(5)) 

H(D)*((D-Vtl) )/(V(2)-V(l) ))**<1.-2.*C<2)> 

A  *( (V(1)*GAM-D)/(V(1)*GAH-V(2) ) )**( 2. *C ( 5 ) -2. *C ( 3) ) 


SBLMAC  3 
SBLMAC  4 
SBLMAC  5 
SBLMAC  6 
SBLMAC  7 
SBLMAC  8 
SBLMAC  9 
SBLNACIO 
SBLMAC 11 
SBLMAC 1 2 
SBLMAC13 
SBLMAC 14 
SBLMAC 15 
SBLNAC16 
SbLMAC 17 
S  BL  MaC 1 8 
SBLMAC19 
SBLMAC  2 C 
SBLMAC  21 
SBLMAC22 
)**C ( SBLMAC  23 
SBLMAC  24 
SDLMAC25 
SBLMAC26 


RSF (D) * AK*( ( CHEN/ ROC ) *D**2) *•=  (1* /  FLOAT (  2+NC  ) ) 
Z1(D)=A*D/(1.^A*D) 

Z«(D)=(B(1)+D (2)*D*A) /( (1.0-D*A)*V( 1)*A) 

ZB(  D)  *AMAX1(-1«  0*  AMINK  1»  0*  ZA(D)) ) 

Z2(D)=EXP(ASIN(ZB(P))/SQRT(B(1)+B<2))) 

>:F(D)-V(1)  +  (V(2)-V(1))*0**(1.i/C(2)) 

IF(MGOOD.tlE.O)  GOTO  10 

MBAD'll 

PRINT  11*NBAD 

RETURN 


11  FORMAT ( lHOi  10X*  30HRETUF.N  FROM  SBLMACH  WITH  NBAD-#  13, 
A26H  AND  WITHOUT  COMPUTATIONS.!*/*  1H  *  10X* 
B34HSUBR0UTII1E  SBLPREP  MUST  BE  CALLED  * 

C39H  BEFORE  OTHER  SBL-ROIITINES  CAM  BE  USED.*/) 

10  I1CAD=0 


Z2Vl*EXP(-l»f 7079 63 268 /SORT (B (l)+8(2) )) 
IF(MR(l).iGE»1.0R.NR ( 2).GE.l )COTO  25 
HBAD=1 

12  PRINT  15*  MBAD 
RETURN 

If  FORMAT ( 1HC*  10X*  30HP.ETURM  FROM  SBLMACH  WITH  NBAD=*I2* 
A39H  DECAUSE  ARGUMENTS  ARE  MOT  WITHIN  RANGE) 

25  IF(RZ.i:Q.O..AMD.TZ.GT.O.)GOTO  35 
:F(RZ.GT.O..AND.TZ.iEQ.C.  )GOTO  45 
MBAD=2 
GOTO  12 


IN  THIS  CASE  THE  INITIAL  MODE  IS  AT  THE  EXPLOSION 
35  TRAT= ( Zl( V( 1) )*Z2V1)/(Z1< V<2) )*Z2(V<2))> 

MUP-I 

IF(TRAT.GT.l.i)NUP— 1 
AMUP=TZ/(Z1(V(1))*Z2V1**NUP) 


SBLMAC27 
SBLMAC  28 
SBLMAC  29 
S8LMAC30 
S8LMAC  31 
SBLMAC  32 
SBLMAC33 
SBLMAC 34 
SBLMAC  35 
SBLMAC  36 
SBLMAC37 
S8LMAC38 
SBLMAC  39 
SBLMAC40 
SBLMAC41 
S8LMAC42 
SBLMAC43 
SBLMAC44 
S EL MaC 4 5 
SBLMAC 46 
SBLMAC47 
SBLMAC48 
SBLMAC  49 
SBLMAC50 
SBLMAC5 1 
SEL  MAC  5  2 
SBLMaC  53 
SBLMAC54 
SBLMAC 5 5 
SBLMAC  56 
SBLMAC  5  7 
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TBUP-TZ 

TFUP*AMUP*Z1( V( 2) )*Z2(V(2 ) ) **NUP 
RFUP*RSF(TFUP) 

AMDOWN*TZ/ { Z1  (V(l )  )*Z2V'1**(-NUP) ) 

TBDOIJN-TZ 

TFD0WN*AI1D0WN*Z1{  V(2)  )*Z2  (V(2)  )**(-NUP) 
RFDOWN-RSF(TFDOUH) 

C  THIS  ESTABLISHED  END  POINTS  OF  BOTH  CHARACTERISTICS. 
C  NOW  GO  TO  FIND  OUT  WHICH  SEGMENT  SHOULD  BE  COMPUTED 
GOTO  55 

95  TS *SQRT( ROC/C HEN ) *(  RZ/AK)**(0.5*FL0AT(2+HC) ) 

C  III  THIS  CASE  A  POINT  ON  THE  SHOCK  IS  SPECIFIED 
TRAT* (Zl( V( 1) )*Z2V1)/ (Z1(V(2) )*Z2(V(2))) 

MUP*1 

IFtTRAT.GT.i.DNUP— 1 
ANUP«TS/(Z1(V(2))*Z2(V(2) )**NUP) 

T8UP-AMUP*Z1(V(1) )*Z2VI**NUP 

TFUP-TS 

RFUP-RZ 

AMD0UN-TS/(ZI(V(2))*Z2(V(2) )**(-NUP)) 
T5D0WN»AMD0UN*Z1(V(1))*Z2V1**(-NUP) 

TFDOWN-TS 

RFDOUM-RZ 

C  AT  55  START  WORK  ON  UPWARD  CHARACTERISTIC 

55  :f(hr(i).ie.o)goto  205 

C  BRANCH  TO  DOWNWARD  CHARACTERISTIC 

IF (RFUP.LT.RNIN.OR.TFUP.LT. THIN )COTO  65 
IF{Q.I»GT.lRMAX«OR.frBUP«GT.  THAXJGOTO  65 
GOTO  ’5 
65  NR( 1) *0 
GOTO  205 

C  NOW  FIND  INTERSECTION  WITH  RMIN 
75  XI-O.l 
Fl»iO.I 
X2-1.I0 
T2-TFUP 
F2*p.rup 
ICNT-0 

P.H«AHAX1(RHIN*0.) 

35  X3«X2-<X2-X1)*( F2-RM) /( F2-FI) 

X3»ANAX1(  0, »  X3) 

IF(X3,LE.0.i)T3«ANUP*Zl(XF(X3)  )*Z2V1**NUP 
ir(X3.GT.0.UT3-AMUP*Zl(XF(X3) )*Z2(XF( X3) ) **NUP 
F3«RSF (T3)*Y ( XF (X3) ) 

IF  (  ABS  ( F3-RM)  •LT«ll»i0E—7*RFD0WN)  GOTO  95 

ICMT-ICNT+l 

:f(icnt.gt.eo)goto  95 

IF(  F3-RM)  86.95*88 
86  X>X3 
F1-F3 
T1-T3 
GOTO  85 
30  X2«X3 
F2-F3 
T2-T3 
GOTO  85 
95  X3*XF ( X3) 


SBLMaC  58 

SBLMAC59 

SSLMAC  60 

53LMaC61 

SBLMAC62 

S3LMAC63 

SBLMAC64 

S3LMAC65 

S3LMAC66 

SBLMAC67 

SBLMAC68 

SBLMAC69 

SBLMAC70 

SBLMAC71 

SbLHAC  72 

SBLMAC73 

SBLMAC74 

SBLMAC75 

SBLMAC76 

SBLMAC77 

SBLMAC 78 

SBLMAC79 

SGLMAC80 

SbLHACei 

SBLMAC82 

SBLMAC83 

SBLMAC84 

SBLMAC 8 5 

SBLMAC86 

SBLMAC 8 7 

SBLMAC88 

SBLMAC89 

SBLMAC90 

SBLMAC91 

SBLMAC92 

SBLMAC93 

SBLMAC94 

SBLMAC95 

SBLMAC96 

SBLMAC97 

SBLMAC98 

SBLMAC99 

SBLMA100 

SBLMA101 

SBLMA102 

SBLMAI03 

SBLMA104 

SBLMA105 

SBLMA106 

SBLMA107 

S6LMA106 

SBLMA109 

S6LMA110 

SBLMA111 

SBLMA112 

SBLMA113 

SBLMA114 
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:F{T3-TMAX}  115/105,65 
105  NP.(1)«1 
GOTO  135 

115  lF(T3*GE<iTMIM)  GOTO  135 

C  BRANCH  IF  FIRST  NODE  UAS  FOUND.  ELSE  GET  INTERSECTION  WITH  THIN 
X1  =  X3 
F1»T3 
>:  2  *  V  C  2 ) 

F2-TFUP 

T2=TFUP 

ICMT-0 

125  X3«X2-(X2-X1) *{ F2-TMIM) / ( F2-F 1) 

T3«=AMUP*Z1<X3)*Z2  <X3)**NUP 
F3-T3 

:F<ABS<F3-TMIN).L7.ii.0E-7*TBUP)  COTO  135 
ICNT-ICNT+1 

ifccnt.ct.zojgoto  135 

IF{F3-TM:N)  126/135/126 

126  X1=X3 
Tl*T3 
F1=F3 
GOTO  125 

128  M2=X3 
T  2=T3 
F2  =  F3 
GOTO  125 
135  VIN=X3 

P.S«RSF  (T3 ) 

T<1,1)«T3 
P.(1/1)*R5*Y(  VIN) 

IFtP.tl/l)  .GT.RHAXJGOTO  65 
USMCKX  2«i/FL0AT  (2+MC ) )  *RS/T  (1/1) 

UP{  1,  U* (R(l,  1) /T(i/ 1 )  }*VIM 
P 1 1/ 1 )  =  1 2.  /  { GAM+1. ) )  *P.0C*USHCK**2*G  ( VIN ) 

F.011/1)  =(  <GAM+1.)/(GAM-1.  ) )  *F.OC*H  ( VIN ) 

DP  (1/1)*R0<1,1}*UP<  1,1)  **2*0.15 
IF(NP.(1)  .EQ.DGOTO  205 

C  BRANCH  TO  COMPUTATION  OF  DOWNWARD  CHARACTERISTIC 
C  MOW  FIND  END  POINT  OF  CURVE 

IF(R(1/1).EQ.RNAX)G0T0  177 

;:><<vin-v<i)  )/<y<2)-v<im**c{2) 

Fl-Rtl/l) 

X2-I.I 

T2=TFUP 

F2=RFUP 

ICMT'O 

RM=AMIN1  ( RHAX/P.FUP) 

145  X3=X2-{X2-X1)*{ F2-RM) /( F2-F1) 

T3* AMUP*Z1( XF ( X  3 ) )*Z2(XF(X3)) **NUP 
F3=RSF(T3)*Y(XF(X3) ) 

ir(ABS(F3-RN) .LT.1,OE-7*RFDOWN)  GOTO  155 
ICNT-ICNT+l 

:f(icnt.gt.20)gdto  155 

I F  C F3-RN)  146/155/148 

146  X1  =  X3 
T1=T3 
Fl=  F3 


SBLMA115 

SBLMAllt 

SBLMA117 

SBLMaIIB 

S6LMA119 

SBLMA120 

SBLMA121 

SBLNA122 

SBLMA1 23 

SBLMA124 

TBLMA125 

SBLMA126 

3ELHA127 

SBLMA12  6 

SBLNA129 

SBLMA130 

SBl.MAl  31 

SBLMA132 

SBLMA133 

SBLMA134 

SBLNA1 35 

SbLMAl 36 

SBLNA1 37 

SBLMA138 

SBLMA1 39 

SBLMA140 

SBLMA141 

SBLNA142 

SBLHA143 

SBLMA144 

SBLMA1A5 

SBLMA146 

SBLMA147 

SBLMA148 

SBLMA149 

SBLMA150 

SBLMA15 1 

SBLMA152 

SBLMA153 

SBLMA154 

SBLMA155 

SBLMaI 56 

SBL  MAx57 

SBLMA158 

S3 L MAI  59 

SBLMA160 

SBL MAI 61 

SBLMA162 

SBLMA163 

SBLMA164 

SBLMA165 

SBLMA166 

SBLMA167 

SBLMA168 

SBLMA169 

TSLMA1 70 

SBLNA171 
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GOTO  145 
148  X2-X3 
T2-T3 
T2«F3 
GOTO  145 
155  X3*XF  (X3) 

IF(T3*LE*TMAX)  GOTO  175 

C.  BRANCH  IF  END  OF  CURVE  FOUND.  ELSE  FIND  INTERSECTION  IJITH  TMAX 
Xl-VIN 
F1*T( 1) 

X2-X3 

F2*T3 

ICIIT-O 

165  X3-X2-(X2-Xl)*(F2-TMAX)/(F2-ri) 

T3«AMUP*Z1(X3)*Z2(X3)**NUP 

F3«T3 

IF  (  ABS  ( F3—TNAX)  * LT.il* OE-7  *THAX)GOTO  175 
:CNT*ICNT+1 

IF( ICNT.G7.20JG0T0  175 
IF{F3-THAX)  166.175/168 

166  ;:i«X3 
T1=T3 
F1-F3 
GOTO  165 

168  X2-X3 
T2=T3 
F2«F3 
GOTO  165 
175  VEN«X3 

:f(veii.gt.v:n)goto  185 

177  HF.{  1) -1 
GOTO  205 
135  KUP-NR(l) 

XIM*( { V I n-V ( 1 ) ) / C V t  2 ) -V C 1 ) ) )**C<2) 

XEM“(  (VEN-V'Cl)  )/(V(2)-V(l))  )**C(2) 

DO  195  KA*2/KUP 

X=XIN+( XEN-XI N) *F  LOAT(KA-l) /FLOAT (KUP-1 ) 

x=xr (X) 

»AHAX1(X/V  ( 1 ) ) 

I F { X.LE.V (1 ) )  T(1/KA) =AMUP*Z1 (V(l ) ) *Z2V1**NUP 
IFtX.GT.Vd) )  T  (1/K.A)  *AMUP*Z1  (X)*Z2{X)**NUP 
RS*RSF(T{ 1/KA) ) 

P.(1/KA)*RS*Y ( X) 

USHCK*  ( 2./FL0AT  { 2+NC )  )*P5/T  (1/K.A) 

P< 1/KA) «  (  2.  /  { GAM+l.l) )  *R0C*USHCK**2*G<  X) 

UPd.  KA)  =  X*P.{1/ KA)/T(1/KA) 

RQ( I. KA) ■ ( ( GAM+1* )/ (GAM-1. ) )*ROC*H( X) 

DP{1/KA)*R0{1/KA)*UP( 1/KA  )**2*0. 5 
195  CONTINUE 

C  AT  205  START  COMPUTATION  OF  DOWNWARD  CURVE 
205  IP<NR<2).LE.0)RETURN 

IF{  RFDOWN.  LT.iRMIN.OR.  TFDOWN. GT* TMAX ) GOTO  215 
IF { 0«« GT.RMAX. OR.TBDOWN.LT, TM IN) GOTO  215 
GOT  }  225 
215  NP.  (  2)  *0 
RETURN 

C  AT  225  FIND  DOWNWARD  INTERSECTION  WITH  RMIM 


SBLNA172 

SBLMA173 

SBLMA17A 

S3LMA175 

SBLMA176 

SBLMA177 

SB L MAI  78 

SBLMA179 

SBLNA180 

SBLHAiei 

SBLMA132 

SBLMA183 

SBLNmI 84 

SBLMA185 

SELMA  166 

SBLMA187 

SBLMA188 

SBLflAi  89 

SBLMA190 

SBLMA191 

SBLHA192 

SBLMA.193 

SDLMA194 

SBLMA195 

SBLHA196 

SBLMaI 97 

SBLMA198 

SBLMA199 

SBLMA200 

SBLMA201 

3BLMA202 

SBLMA203 

SBLMA204 

SBL.MA205 

SBLMA20O 

SBLMA207 

SBLNA208 

SBLMA209 

SBLMA21C 

SBLMA21 1 

SBLMA2 12 

SBLMA213 

SBLMA214 

SBLMA215 

SBLMA216 

SBLMA217 

38LMA218 

SBLMA219 

SBLMA220 

SBLMA221 

SBLMA222 

SBLMA223 

SBLMA22  4 

SBLMA225 

SBLMA226 

SBLMA227 

SBLHA228 
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225 

x>o.o 

SBLNA229 

ri-o 

SDLMA230 

X2-1.0 

5BIMA231 

T2-TFD0UII 

SBLNA232 

F2=  RFDOWU 

SbLMA233 

ICNT-0 

SDLNA23A 

F;M»AM«Xl(RMIN/0«  ) 

SBLMA235 

235 

X3=X2-(X2*-X1 ) *{ F2— RM J /( F2—F1) 

SBLMA236 

::3*AN/,X1(0.»X3) 

SDLMA2  37 

IF ( X3»lLE«0«i)  T3»AHD0WN*ZI(  XF  ( X 3 ) ) /Z2V1**NUP 

SBLMA238 

IF  (  X3»CT»0*i)  T3»AND0WM*Z1(  XF(X3))/Z2(XF(X3))  **NUP 

SBLKA239 

F3*P.CF  (T3)*Y (XF  (>'3J ) 

SBLMA2  40 

IF(ABS(F3-F:tl)  .LT.il.iOE-7  *RFDOWM) GOTO  245 

SBLMA24i 

ICHT-ICHT+1 

S6LMA242 

IF  { ICNT«GT*i20  )GOTO  245 

SBLHA243 

:F(F3-RHJ  236/245/230 

SELNm244 

236 

X1  =  X3 

SBINA2A5 

T  I=T3 

SBLNA246 

F1=F3 

SBLMA247 

GOTO  235 

S3LNA248 

238 

.-I2-X3 

SBLMA249 

72*73 

SBLMA2  50 

F2=  F3 

SbLMA251 

GOTO  235 

SBINA252 

24? 

X3*Xr(X3) 

SBLMA253 

IF ( T3-TNIN) 215/ 255/ 265 

SBLMA254 

2  55 

NR( 2) =1 

S  B  L  HA 2  5  5 

GOTO  285 

SBLMA256 

265 

IF(T3«LE#TMAX)G070  285 

SBLI1A257 

BRANCH  IF  FIRST  NODE  FOUND 

S8LNA258 

X1=X3 

S6LMA259 

F1  =  T3 

SBLNA260 

X2=V( 2) 

SBLMA261 

T2-TFD0WN 

SBLMA26  2 

F2=T2 

SBLMA2  63 

ICNT=0 

SBLMA264 

275 

X3=X2-(X2-X1) *( F2-TMAX) /( F2-F i) 

SBLHA265 

T3=AMD01JII*Z1(X3)/Z2(X3)**NUP 

SBLMA266 

F3  =  T3 

SELHA267 

IF(  ABS  (F3-TNAX)  »  LT«I1«  OE-7  *TBD0WN)G0T0  285 

SBLHA2  68 

ICNT-ICNT+1 

SBLMA269 

I F { ICNT»GT«20)GOTO  285 

SBLMA270 

:F(F3-THAX)  276/285/278 

SBLMA271 

276 

X2=X3 

SBLMA272 

T2=T3 

SBLMA273 

F2=  F3 

SBLMA274 

GOTO  275 

SBLMA275 

278 

Xi*X3 

S6LNA276 

T1*T3 

SBLHA277 

Fl=  F3 

SBLMA278 

GOTO  2^5 

SBLMA2  79 

23? 

V  IN*X3 

SBLMA280 

FiS-RSF(T3) 

5BLNA281 

T(2/1)=T3 

SBIHA282 

r.(2/i)=r.s*Y(viN) 

SBLMA2  83 

IF(R(2/1).GT.RMAX)G0T0  215 

SBLMA284 

USHCK*(2«/FL0AT(2+NC)  )*RS/T(2/1) 

SBLHA2e5 

63 


P  (2.  1 ) « ( 2*1/  ( CAM+1* ) )  *R0C*USHCK**2*G  (VIN) 
UF(2/1)«VIN*R(2/1)/T(2/1) 

R0(2/ 1) *(  (GAN+l*)/( GAM-1*  ) )*ROC*H(VIN) 

DP  ( 2/ 1)  «R0(  2/  1)  *UP(  2/ 1)  **2*0*15 
IF(R(2/1),EQ*RMAX)NR(2)«1 
IF(MR(  2)  .IEQ.1  )RETUP.M 
C  :!0U  FIND  END  POINT  OF  DOWNWARD  CURVE 
C  FIRST  FIND  INTERSECTION  WITH  P.MAX 

XI* ( ( VIN-V( 1 ) )/ (V(2)-V(l) ) )**C(2) 
n«R(2/l) 

X2-1.0 

T2-TFD0WN 

F2-RFD0WM 

ICNT-0 

RM*  AI1INK  RFDOWM/RMAX ) 

295  X3*X2-(X2-X1) *( F2-RM) / ( F2-F1) 
T3«AMD0'JI1*Z1(XF(X3))/Z2(XF(X3))**NUP 
F3«RSF(T3)*Y(XF(X3)) 

IF (  ABS  ( F3— RM )  * LT#ll«OE— 7  *RFDOWN)GOTO  305 
ICNT-ICMT+1 

IF(ICMT.GT.20)C0T0  305 
I F  (  F3-P.M)  296/305*298 

296  X1*X3 
T1  =  T3 
Fi«r3 
GOTO  295 

29C  X2«X3 
T  2=73 
F2-F3 
GOTO  295 
305  X3-XF(X3) 

IF (T3«GE«TMIN )GOTO  325 

C  BRANCH  IF  END  POINT  FOUND.  ELSE  FIND  INTERSECTION  WITH  THIN 
XI«VIN 
Fi-T(2i  1) 

F2=T3 
X2  =  X3 
ICNT*0 

315  X3=X2-(X2“X1 ) *( F2-TMIN) / ( F2— FI) 

T3«AMD01M*Z1 ( X3 )/ Z2 ( X3) **NUP 
F3*T3 

IF  ( ABS  ( F3-TIIIN) »  LT*il*0L— 7  *TBDOWM)GOTO  325 
ICNT*ICHT+1 

IF ( :CNT.GT*2O)G0T0  325 
IF  (  F3— TNIil)  316.  3 25/  3 IB 

316  X2*X3 
T2=T3 
F2-F3 
GOTO  315 

31C  X1*X3 
T1=T3 
F1  =  F3 
GOTO  315 
325  V2U-X3 

IF ( VEN*GT«VIN)GQTO  335 

NR( 2) *1 

F.ETUF.N 


S&L.MA286 

SOI HA287 

S3LMA2db 

5BLMA289 

S&LMA290 

SBLMA291 

SBLNA292 

S8LHA293 

SbLH«294 

SBLMA295 

SCLHA296 

SDU1A297 

SBLMA298 

S3LNA299 

SBLMA200 

SBLNA3G1 

SBLMA302 

SBLHA303 

SB  L  nA  304 

SBLMA3G5 

SBLMA306 

SBLHA307 

SBLMA308 

SBLNA309 

SLLHA310 

SBLMA31 i 

SBLHA312 

SELHA313 

SGLNA314 

SBLHA315 

S BLN A 31 6 

S6LMA317 

S6LHA318 

S  B  L  NA  319 

SBLMA320 

S&LMA321 

SLLHA322 

SBLNA323 

SBLHA324 

SBLHA325 

SBLHA326 

SBLMA327 

SBLHA328 

SBLHA329 

SBLHA33G 

S3LMA331 

SBLHA332 

S3LHA333 

SBLHA334 

SBLHA335 

SbLHA336 

SDLMA337 

SBLMA333 

SBLMA339 

SBLMA34C 

SBLMA341 

SBLHA342 
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C  LOOP  TO  COMPUTE  DOWNWARD  CHARACTERISTIC 
33?  KUP-MR< 2) 

XIN«( (VIM-V(l) )/(V(2)-V(l)) )**C (2) 

X2N“(  ( VcM-V (1))/(V(2)— V'(l))  )**C ( 2 ) 

DO  3 A?  KA-2/KUP 

X=XIM+ ( XEM-XI M) *FLOAT(KA-l) /F  LOAT (KUP-1) 

x-xrm 

X-ANAX1(X/V(1)) 

IF(X.LE.V(1) )  T(2/KA)*=AMD0WN*Z1(V{1)  )/Z2Vl**NUP 
:r  (X.GT.V(l) }  T(2/KA)  “AMDOIJN* Z1  ( X. ) / Z2 { X ) **MUP 
P.S«RSF(T(2/KA)) 

R ( 2/KA) =RS*Y (X) 

U:HCK«(2.I/FL0AT(2+NC)  )*P.S/T(2/KA) 

P  (2/KA) * (  2*1/  (  GAM+1«I) ) *ROC*U3HCK**2*G ( X ) 

UP(  2/ KA)  *X*P.(  2/ KA) /T (2/KA) 

PO( 2/ KA)  = ( (GAM+1, )/ (GAh-l.) )*ROC*H(X) 

DP( 2/ KA) -RO{ 2/KA) *UP ( 2/KA  )**2*0.5 
3 A?  CONTINUE 
RETURN 
END 


SBLMA3A3 
SBLMA344 
SBLf A3A5 
SBLMA346 
SBLNA3A7 
SBLMA348 
SbLNA3A9 
SBl MA35C 
SBLMA351 
SBLMA352 
SBLMA353 
SBLNA354 
S £j l  HA35? 
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USER  EVALUATION  OF  REPORT 


Please  take  a  few  minutes  to  answer  the  questions  below;  tear  out 
this  sheet  and  return  it  to  Director,  US  Army  Ballistic  Research 
Laboratory,  ARRADCOM,  ATTN:  DRDAR-TSB,  Aberdeen  Proving  Ground, 
Maryland  21005.  Your  comments  will  provide  us  with  information 
for  improving  future  reports. 

1 .  BRL  Report  Number 

2.  Does  this  report  satisfy  a  need?  (Comment  on  purpose,  related 
project,  or  other  area  of  interest  for  which  report  will  be  used.) 


3.  How,  specifically,  is  the  report  being  used?  (Information 
source,  design  data  or  procedure,  management  procedure,  source  of 
ideas,  etc.) _ 


4.  Has  the  information  in  this  report  led  to  any  quantitative 
savings  as  far  as  man-hours/contract  dollars  saved,  operating  costs 
avoided,  efficiencies  achieved,  etc.?  If  so,  please  elaborate. 


5.  General  Comments  (Indicate  what  you  think  should  be  changed  to 
make  this  report  and  future  reports  of  this  type  more  responsive 
to  your  needs,  more  usable,  improve  readability,  etc.) 


6.  If  you  would  like  to  be  contacted  by  the  personnel  who  prepared 
this  report  to  raise  specific  questions  or  discuss  the  topic, 
please  fill  in  the  following  information. 

Name : 


Telephone  Number: 
Organization  Address: 


